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SECTION  1 


RADIO  WAVE  PROPAGATION  IN  STRUCTURED  IONIZATION 
FOR  SATELLITE  AND  RADAR  APPLICATIONS 


1.1  INTRODUCTION. 

This  report  is  an  extension  of  References  and  2*  which  developed  the 
algorithms  necessary  to  calculate  the  signal  structure  parameters  and  the  generalized 
power  spectrum  which  represent  the  simultaneous  phase  and  amplitude  effects  of  prop¬ 
agating  electromagnetic  waves  through  media  characterized  by  a  structured  index  of 
refraction.  Those  earlier  works  were  restricted  to  single  component  power  spectrums, 
did  not  treat  antenna  effects  with  sufficient  generality,  and  did  not  treat  phase  oply 
(total  electron  content)  effects  adequately  for  space  radar  applications.  This  paper 
will  eliminate  these  restrictions. 


1.2  ENMRONMENT  CHARACTERIZATION. 

The  first  sten  in  any  propagation  study  is  to  characterize  the  ionization 
)  i  .-Mvalsntiy  the  index  of  refraction  fluctuations  of  the  ionospheric  propagation 
environment.  Figure  1-1  illustrates  the  geometry  of  a  typical  satellite  link.  The  index 
of  refraction  fluctuations,  represented  schematically  by  the  curved  lines  are  typically 
elongated  along  the  magnetic  field.  The  unit  vector  t  is  parallel  to  the  magnetic  field 
is  a  slowly  varying  function  of  space  since  the  field  lines  are  curved.  The  i  axis  and 
the  i  unit  vector  represent  the  propagation  line  of  sight.  One  end  of  the  radio  link  is 
at  2  =  0  and  the  other  end  is  at  2  =  Except  where  specifically  noted,  all  quantities 
calculated  are  at  2  =  2|.  The  f  and  s  unit  vectors  complete,  with  f,  an  orthogonal 
coordinate  system  that  is  used  to  define  the  structure.  The  orientation  of  f  and  I 
is  chc^en  to  best  represent  any  anisotropy  of  the  index  of  refraction  structure  about 
the  field  line.  As  with  f  and  s  may  be  slow  functions  of  position.  In  the  r,  I,  and 
i  system,  the  index  refraction  fluctuations  are  represented  by  a  power  law  power 
spectral  density. 

‘Wittwer,  L.  A.,  Redip  Wave  Propagaitoa  tii  SlruefareJ  hniiatioa  for  SaUlliU  Applications,  DNA 
5304D,  31  December  1979. 

^Wittwer,  L.  A.,  Radio  Propagation  Structured  lonisaUon  for  Satellite  Applications  ll, 

DNA-IR-«3^2,  1  August  1982. 
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Figure  1-1.  Pr(^>agatloD  enviroument. 
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tn  ttt 


V^nlK  K  V  <  8.r»/»A»?W,(n,  n'.B)L,L.J,.r(n)/r(n  -  3/2) 

(1.1) 

where 


Xr>  X«,  Xt  =  structure  outer  scales 

ifit-tM  =  structure  freezing  scales 

An|  =  index  of  refraction  variance 

r(n)  =  gamma  function  of  argument  n 

2n  -  2  =  intermediate  scale  spectral  index  (1.5  <  n  <  2) 

2a'  -  2  =  transition  scale  spatial  index  (2  <  a'  <  4) 

It  is  assumed  that 


X,  X. 


^-R<1 

iH 


Th^  assumptions  regarding  the  f,  s  and  £  axes  greatly  simplify  the  following  de¬ 
velopment  without  signihcantiy  limiting  the  applicability  of  Equation  1.1.  It  is  not 
claimed  that  Equation  1.1  with  R  invariant  fully  reprints  all  aspects  of  plasma  or 
other  structure.  Indeed,  it  does  not.  It  does,  however,  provide  an  adequate  represen¬ 
tation  for  th(»e  specific  features  that  impact  the  distortion  of  radio  or  radar  signals. 
Equation  1.1  has  also  found  application  in  specifying  infrared  and  other  clutter  back¬ 
grounds.  With  the  above  assumption 


N,(n.n'.fi)  =  fl  -  ( 

where  h{^n\R)  b  defined  in  Appendix  £  in  Subroutine  PROP. 

The  structure  variation  of  the  index  of  refraction  perpendicular  to  the  i  axis 
dominates  the  propagation  effects  while  the  variation  parallel  to  the  s  axis  enters  only 
through  the  strength  of  the  integrated  phase  variance.  Thus  Equation  1.1  must  be 
transformed  to  a  frame  with  one  axis  being  Utc  B  axb-  These  new  axes  arc  defined 
by: 
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(1.2a) 


ixz 

\ixB\ 


zx{ix  z) 
ix{tx  5)1 


(1.2b) 


Two  rotations  suffice  to  accomplish  the  transformation.  First,  the  f  axis  is  rotated 
about  the  t  aids  by  angle  ^  to  become  parallel  to  the  5  axis.  The  angle  4>  Is  defined 
by 


isia*  = 

It  X  5| 

(1.3a) 

^  r .  (t  X  5) 

cos4>  =  }  ' 

\txz\ 

(1.3b) 

Next,  the  t  axis  b  rotated  about  the  z  axis  by  angle  ^  into  the  z  axis.  The  angle  4 
b  defined  by 


xsin^  fx5 

cos#  =  i-i  (l*4b) 

These  transformatimis  can  be  simplified  by  defining  new  effective  scale  sizes. 

Ll  =  L*  cos’ #  + £f*  sin*  #  (1.5a) 

L*  -  (Lj  sin* #  +  Lj cos* cos*  #  +  sia*  4*  (l*5b) 

Ll  =  [LI  sin*  #  -f  L*  cos*  #)  sin*  ♦  -f  Lj  cos*  #  (1.5c) 

=  (L*  -  L*)cos#co8^5in#  (l.5d) 

s:  (i»*  -  Lj)sin#cosdsin#  (l-Se) 

L^,  =  m  -  L]  sin*  d  -  cos*  ^)  c(»  #  sin  #  (1 .5f) 


Simitar  transformations  apply  to  the  freezing  and  inner  scales.  The  final  result  in  the 
X,  y,  z  coordinate  system  b 
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PSD(K„  Ky,K,)  =  SK^f^An^Ni{n,n\R)LrL,LtT{n)/T{n  -  3/2)  / 

[(1  +  Kill  +  Kill  +  KlLl  +  2L^yK,Ky  +  2L„K^K,  +  2Ly,KyK,)^ 

(1  +  KUI  +  Klti  +  Kill  +  2l,^K^Ky  +  2t^K^K,  +  2ty,KyK,f-*]  (1.6) 

where  the  small  t  parameters  are  the  transformed  freezing  scales  The  scales  and 
are  not  signal  duxorrelation  distance  in  this  report. 


1.8  DERIVATION  OF  SIGNAL  STRUCTURE  PARAMETERS. 


The  required  signal  parameters  are  derived  from  the  differential  phase  spec¬ 
trum  which  is  Equation  1.6  with  =  0  and  multiplied  by  K*  where  K  is  the  carrier 
frequency  wave  number.  See  Appendix  D  in  DNA  S304D.  The  differential  phase  sp^- 
trum  is 


dP^{K^Ky)  ^  > An?iV5(n,  n\  R)LrL,LtT{n)/T{n  -  3/2)  / 

dz 

[(1  +  Kill  +  Kji;  +  2L^K.K,)'‘{1  +  Kill  +  KX  + 

(1.7) 

The  mean  square  phase  fluctuation  b 


where 


daj 

dz 


(1.8) 


dcrl  2y/^'(n  -  l)Ns(nyn\R)K^  An*  L' 
dz  ~  r(n-3/2)/V,(n,n',R) 

jVi(n,B',R)  =  [■  “  (^rry)  (/>(»•’»'•**)*)’"’ 


(1.9) 

(1.10) 
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where  /j(n,n',i2)  is  defined  in  Appendix  E  in  Subroutine  PitOP.  The  scale  is  the 
effective  z  axis  scale.  Finally 


=  §  4»^L>i5  -  Li,{n  -  l)N,(n,  n',  R)  / 

[(l  +  Kill  +  KjLj  +  2L.,Jf.if,)”  (l  +  Kie,  +  Kit?,  +  2<,.Jr.K,)"'‘"] 

(1.11) 

In  many  applications  it  is  necessary  to  know  the  one  dimensional  power  spectral 
density  along  a  fixed  direction  in  the  x,  y  system.  Let  the  direction  be  defined  by  the 
vector,  xq  i  +  yoj.  The  spectrum  is  approximated  by 


dP^jK)  _dcl _ 2^N2{n^',R)T{n  -  1/2){1  +  6.4/n')  L 

dz 


where 


iz  r(n  -  1)(1  +  6.4/»)  (1  +  (d?  + 


(1.12) 


=  R^^  +  {1-R^^)gkp 


-6.4 


Ln  -t-6.4nj 
{xl  +  yimLl-LD 
^IVo  +  ~  2L*v^oyo 

£2  =  R^L^=. 


n'  7^  n 


n  —n 


4i!/o  +  <5*0  -  2«„Io!A) 
Let  me  now  define  the  “local”  phase  correlation  function: 


1 

^^(p(a:o,yo))  =  5^  /  d^cos(ifLp(xo,yo)) 

J-oo 


dz 


1^ 

‘  dz 


where 


p^(xo,yo} 


^xVo  "I"  "  2LjyXot/-) 


xp  +  yg 

L» 
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Finally,  let  me  define  the  “local”  phase  structure  function: 


=  l-i2^(p(xo,yo)) 


(1.13a) 


27r  J-o 


dK 


2  sin 


,  (KLp{io,y^)\  dPt(K) 
\  2  )  iz 


dal 

j 

’  dz 


(1.13b) 


Appendix  C  documents  a  Fortran  function  routine  to  evaluate  the  phase  structure 
function. 


We  now  consider  the  calculation  of  signal  parameters  with  reference  to  a 
fixed  coordinate  system  perpendicular  to  the  z  axis  The  local  x,  y  system  can  rotate 
in  some  situations  as  a  function  of  z.  Let  us  define  the  u  and  v  axis  perpendicular  to 
the  z.  At  each  point  on  I'..''  B  axis,  the  rotation  angle  between  the  x  and  u  axes  is 
defined  by 


z  sind  =  St  XU 

COS$  =  X'U 

In  terms  of  the  system 

,  _  Llul  +  Llvl  -  2L^uoVo 
P  (tio.%  - 

"uv 

where 

L*  =  L*  cos*  6  LI  sin*  9  +  cos  6  sin  9 
L\  =  L\  sin*  9  ■¥  L\  cos*  9  -  2L,^  cos  9  sin  9 
-  {^l~  ^l)  sin  5  cos  d  +  Lsj,  (c<»*  9  -  sin*  5) 
The  local  uq,  vq  coordinates  along  the  line  of  sight  are 


=  u(0  -s-vw. 
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where 


V  =  -  ^  i^Los(O)  +  ^ST  -  ^LOs(«t) 

t4os(«)  =  line  of  sight  velocity 
Vsr{z)  =  structure  velocity 

The  primed  coordinates  are  at  s  —  0  and  the  unprImed  coordinates  are  ht  z  —  Zt. 
The  complex  signal  space-time  correlation  function  is 

R{u,v,u',v\t)  =exp|- 

For  most  applications  this  function  will  have  the  form 


i2(u,v,tt',v')  =  exp 


I- 


rtt*  t*  tiv  ut  v< 

^  ^  r|  “  “"tx  ~  “Ob  "  “Ob 


u'®  w*’ 


u'v' 


ttu'  t»v' 

—  2Cm,»  j-T  2C7ut,i .  . 


-  similar  cross  terms 


(1.15) 


where  only  the  most  important  terms  have  been  written  explicitly.  For  convenience, 
let 

I(u,v,u\v\t)  =  ^  d3~D^{p(u<,(u,u',(),Uo(u,u'.0)l 
Now  4i.  Ai>  and  i^,t  and  Tq  are  found  by  solving 


;(£u.0, 0,0,0)  =  C(al) 
/(0.4, 0,0,0)  -  C(al) 
/(0,0,4s0.0)  =  C{a5) 
7(0, 0,0,  £,,,0)  =  C(aJ) 
7(0,0, 0,0,ro)  =  C{o*) 
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The  cross  terras  require  a  different  procedure.  For  exaraple,  consider  Cut,. 

•  Solve  0,0,0)  =  C(oJ). 

•  Let  =  1  -  C(<7})/(2o*). 

•  If  Cuv  >  —0.26,  then  you  are  finished. 

•  Otherwhe,  solve  /(-Oi^,aiA,»0,0,0)  =  C((7j). 

•  Then  =  -1  +  C(o})/(2a»). 

where 

'  (i  - 

For  defining  “reasonable  worst  case”  pa-araeiers,  we  seek  the  minimum  decorrelation 
lengths  at  a  =  0  and  a  =  a*.  At  “  =  define  e 


U  £u  >  ivt  then  let  c  =  €  +  5(i/2.This  angle  defines  a  rotated  coordinate  system  with 

u  X  p  =  a  sin  f 
u  p  -  cos  t 

p  =  u  cos  f  +  V  sin  £ 
g  —  -u  sin  £ -f  V  cost 

The  new  correlation  quantities  of  dominant  concv;rn  are 


i  2  ,  *  •  2 

—  cos  £  -  — —  COS  V  Sin  £  +  sm  £ 


{1.16a) 
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(1.16b) 


1  1  .  ,  2C„«  .  1  . 

^  =  ^sm*t  +  -^co8C8ine+^cos*£ 


C'p. 

=  0 

Cpi 

=  41 

c,t 

=  4| 

1  « 

A, 

Cti* 


The  energy  angle-of-arrival  variances  are 


(1.16c) 

(1.16<i) 

(1.16e) 


2 

KHl 


2 

KH\ 


will  always  be  greater  than  or  equal  to  af^,  so  it  is  the  “reasonable  worst  case" 
value.  The  probabiUty  density  function  for  arrival  angles  is  proportional  to 


p(s„e,) 


or  in  “wave  number"  space 


1 


exp 


P{,K„K,)  =  ^(,1,  exp  [-i  (k;^  +  Kit]) 


(1.17) 


The  same  general  set  of  procedures  for  calculating  “reasonable  worst  case"  parameters 
must  be  executed  at  a  =  0  since  the  spatial  parameters  are  not,  in  general,  reciprocal. 
At  the  transmitter  (a  =  0),  we  are  generally  content  with  calculating  ^d  (or 
and  <7|^). 

At  the  receiver,  the  signal  correlation  function  is  written  as 


The  corresponding  power  spectrum  is 


F{K,,K„J) 


JVo 


-Cj) +  «•,’/;(! -Ci)+ 


(2jrro/)’  +  2ifpVf,^Cp,C,«  -  "  2if,£,(2;rro/)C,t]} 


(1.18) 


where  Nq  =  ^1  -  C’j  -  C*|  . 

For  radar  application,  it  is  necessary  to  know  the  decorrelation  between  the 
outward  path  to  the  target  and  the  return  path.  Assuming  that  the  radar  is  at  2  = 
the  two  way  decorrelation  is 


I  (1-19) 

where 

c  =  speed  of  light  =  3  x  10*m/sec 
Vrx  = 

Vj>y  =  Vt(2)  •  y 

VVW  =  VstW  - 

We  define  the  parallel  signal  decorrelation  time  as 


TWPD  =  exp 


{-/o 


^  ,  dal 
dz  Da 
dz 


2z 

c 


\ 


lFt’x  + 


li  -  11 

**^*'*'V  "xy 


_  _ _  2x 

"  X  [  |(vUo  •  2)  -  V.la?  +  •  j)  -  v;^o]  ] 


dal  , 


2) 


The  final  decorrelation  time  is 
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(1.20) 


ib  =  min(ro,ni) 

where  the  Tq  within  the  parentheses  was  calculated  earlier. 

The  next  quantity  required  is  the  frequency  selective  bandwidth,  /q. 


where  Ci  =  delay  parameter  («0.25)  and 


Appendbc  D  discusses  the  derivation  of  the  above  equation.  Also 


(1.21) 


=  Lj  cos*  6  +  sin*  6  +  2L*y  sin  6  cos  $ 

=  L*  sin*0 -I- Zrjcos*^  -  2Ljysinficostf 
=  (L*  -  L* )  sin  d  cos  (cos*  0  -  sin*  0) 
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p 


\jl\ltz[zt  -  z) 

ZiL 

DM 

The  final  generalized  power  spectrum  at  the  receiver  b  usually  expressed  as 


=2i^foF{K,,K,J)6 


2irfoT  - 


(Ki  +  K^) 
v'5+q  4 


(1.22) 


where 


f  dxS{x-a)f{x)  =  /(a) 

J^OO 

/oo 

dx  S{x  -  o)  =  1 

"OO 


The  mean  square  log  amplitude  fluctuation  is  used  to  denote  the  presence 
of  amplitude  fluctuations.  It  is 


X2(.)  =  ^.(»,n'.«)(n-l)£  (l+„)"(.  +  «>v) 


(1.23a) 


«  .5(1  -  cos{May)MMi>y)] 

- - :i.'Tn>-n  (l*23b) 


where 


-  2)  Ll  +  Ll 
~Tzr~  ■  2{LILI  - 


z{zt  -  Z)  s/m  -  Liy  +  ill 

Kz,  ‘  2{LlLl~Ll^) 

Appendix  A  discusses  the  development  of  Equation  1.23  and  presents  an  algorithm  to 
evaluate  Equation  1.23b. 
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The  X*  algorithm  is  accurate  to  better  than  ten  percent.  The  quantity  x* 
should  henceforth  be  carried  as  a  primary  quantity  to  measure  amplitude  fluctuations. 
This  function  was  previously  done  by  the  Rayleigh  phase  fluctuation  parameter. 


1.4  EFFECTIVE  SCALE  SIZES  AND  INDICES. 


For  most  applications,  the  eflective  scales  and  indices  are  used  to  calculate 
the  statistics  of  the  environment  produced  Doppler,  Doppler  rate,  group  delay,  etc. 
At  each  position  along  the  line  of  sight,  we  deflne  the  scales  along  the  local  velocity 
vector.  Then 


ju.  ^  {Vl{^)+VyH^))LlLl 
^  ^  Vi{z)Ll  +  V^^{z)Ll 

■  y^wi^y}m 

= 


G(7n)  = 


|vp 

m 

u 

m  +  1  ‘  ' 

2n,  R)  +  F{2n!  -  m  - 

R' 

=  fi/t 

ti  =  inner  scale 

F(n,R) 

= 

n  ^  0 

=  HUR) 

n  =  0 

=  \A?W  +  TO 

fo  =  /  dz 
Jo 


^  ..  <^1  C(2) 

dz  G{0) 


Finally 


A  ^  JL 

IqJq  ^  dz  G(0)  l^l 


(1.24a) 
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1,8  = 

fl  <1 

[/,  h  lU  G(0)  ||7|. 

V,8 

(1.24b) 

Aeff  ~ 

fi 

[loh  dz  G(0)  |K| 

V,ff 

(1.24c} 

n,s  = 

IqJo  dz  G(0) 

(1.24d) 

<ff  = 

IqJq  ^  dz  G(0) 

(l.24e) 

At  this  point, 


must  be  redefined  to  reflect  the  “Doppler  power.”  Let 


R  =  ^ff/  Leir 

R’  - 

n  = 

=  <» 

1^1  =  Va=e/ro 
L  =  L,ff 

in  the  <*(m)  function.  Then 


(1.25) 


For  most  applications,  V,9  can  be  approximated  by 


1.5  SUMMARY. 


The  preceeding  sections  extended  the  results  of  References  1  and  2  to  permit 
the  prediction  of  the  effects  of  structured  index  of  refraction  on  electromagnetic  waves 
to  two  component  power  spectrums.  Six  new  parameters,  La,  and  Aeff)  n«fr  &nd 
n'a,  were  introduced  to  describe  the  total  electron  content  power  spectrum.  One  new 
parameter  was  introduced  for  some  radar  applications.  Finally,  the  Rayleigh  phase 
variance  has  been  eliminated  in  favor  of  the  mean  square  log  amplitude  fluctuation  for 
indicating  amplitude  fluctuations.  Appendix  E  documents  a  FORTRAN  subroutine 
which  calculates  the  described  signal  parameters. 
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APPENDIX  A 


CALCULATION  OF  THE  MEAN  SQUARE  LOG  AMPLITUDE 

FLUCTUATION 


The  equation  for  the  mean  square  log  amplitude  fluctuation  is 


zi-  [‘•j.  r  ^  r 


IKzt 


dP^{K,,Ky) 

dz 


(A.1) 


where  jg  found  in  Equation  1.11.  Equation  A.l  can  be  transformed  (as 

shown  by  Scott  Frasier)  into 


^  /V  r  ■  cos(Af«(2)y))Jo(M»{z)y))  . 


where 


M*(^) 


3(Z(  —  S)  ^<1  + 

.-(j,  - .-)  \/m  -  ij)’ + *1-1, 

Kz,  2(L\L\  -  LI) 


To  avoid  the  Bessel  function  aiid  the  cosine,  the  numerator  in  the  inner  inU^ral  is 
approximated  by 


.S(l-cos{M,{2)y))Jo(Mi(2)y)! 


a* 


{  0.5 


a  <  >/Cj 

y/c[  <  a  <  3.1 
3.1  <  a 


where 


A-1 


a  =  M{z)y 


Mi,) ... 

b  =  ln(0.5/ei)/ln(3.1/VcD 

^  (l  +  fl>/2)» 

‘  1.2  +  4.812, 

Ri  =  M}/Ml 

Equation  A.2  using  the  above  approximation  tan  be  rapidly  integrated  since  the  inte¬ 
grand  is  power  law  over  long  ranges  of  y.  Points  can  be  chosen  to  closely  approximate 
the  integrand  as  a  series  of  power  law  segments.  Each  segment  can  be  analytically 
integrated.  Appendix  E  shows  a  line  of  sight  integrator  that  uses  thb  method. 
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APPEJNDDC  B 

LARGE  TARGET,  LARGE  APERTURE  SIGNAL  STRUCTURES 


The  signal  at  coordinates  received  from  coordinates  {hi,<f>)  can  be 


written  as 


exp  (iKiyJ +  \gi  -  hi\A 

E{z,guht,Ki)  =  U{z,guhuKi) - ^  ^  (B.l) 

V2*+  tsi  -^ll’ 


where 


91  =  xi »  +  yi  J 

hi  =  ttit  +  Vij 


Applying  the  parabolic  approximation 


The  autocorrelation  function  can  be  written  as 


exp{»7f|2)  exp 


iKilgi  -  A, 


(B.2) 


IiK  -  - 

—  [//'fc  +  C-  g-  fl-H-C-fc] 


{HH  +  GG)  ^(kk  +  gg)  -  -  kg 

s  2  8  4 


(B.3) 


G{z,  I  d,  h,  H,  Aik,  K)  =  U*{z,  gu  K  Ki)U{z,  h.  ^2,  K^) 


ffi + h 


A  A 

X;  rj 


xii-x2  V  _y^^yi 

~  2  '  ~  2 


—t 

Q 


92-91  =  xi  +  yj 


X  =  X2  -  li  ,  y  =  92-91 


H  =  h±h  =  ui  T  V] 


Ui  +  U2  y  _  Vl  +  Vj 

2  ’  ”2 

A  ^ 

h  =  h2-  hi  =  ui  +  vj 


tt  —  U2  -  Ui  ,  V  =  Uj  -  Vi 


Ak  =  Xj-Xi 


K  = 


/ri  +  X; 
2 


The  equation  for  G{z,9,G,h,H,Ak,K)  is 


dz 


iAk 


2KxK2 


d^G  d^G  1  p’G  d^G 
dx^  i  vax^  ''■  ay* 


)i 


+ 


iK 

KxK, 


■  d^G  d^G 
dxdX  dydY 


{xj-j^^  (X  -  U)  dG  _  (y  -  v)  dG  _  {Y  -  V)  dG 
z  dx  z  bX  z  by  z  bY 


1 

2 


(B.1) 


where  is  evaluated  for  a  wavenumbex  of  \jK\Ki.  Also,  V*  and  Vy  have  been 
introduced  to  represent  the  motion  of  the  environment  with  respect  to  the  line  of 
sight.  The  mixing  of  the  notation  for  the  wavenumbers  was  for  convenience.  is  the 
normalized  integrated  phase  autocorrelation  function. 


/2^(0,0,0,0,z)  =  1 

Vs  and  Vy  are  local  velocities  and  t  is  the  time  displacement.  The  exponential  term 
divided  by  a’  in  Equation  B.3  contains  the  free  space  solution  for  the  total  signal 
autocorrelation  function.  Thus  in  the  absence  of  any  signal  scattering,  G  is  always 
equal  to  one.  In  general,  does  not  depend  on  X  or  Y.  Since  G  does  not  have  any 
dependence  on  X  or  y  in  free  space,  it  can  never  develop  any.  Equation  B.4  caii  thus 
be  written 


dG{z,g,h,Ak,t)  iAk 


dz 


2K,Xj 


^  dy* 


(z  -  u)  dG  {y  -  v)  dG 
z  dx  z  dy 


Kti-Kl 


2  [  KxKt 


dal  dcr\ 

-^G  +  -  V,t,y-  V,t,z)G 

az  az 


(B.5) 


Let 


£^#(p(^»yo))  ==  1  -  /4(xo,Sto,2) 


B-3 


where 


+  Llvl  -  2L,yXoyo 

LlLl  - 

environment  outer  scales 

The  z  dependence  in  is  implicit  in  the  scale  sizes  and  the  exact  functional  depen¬ 
dence  of  R^.  is  the  “local*  phase  structure  function. 

Equation  B.5  simplifies  to 


p\xo,yo]  = 

^iit  ~ 


dG{z,lh,Ak,K) 

dz 


iAk  \d^G  d^G 
2KxKi  [ai»  ay* 


(g  —  u)  dG  (y  -  v)  dG 
z  dz  z  dy 


2K1K2  dz 


do^ 

dz 


DM^-y^Uy-yy^))G 


(B.6) 


Let 


G  —  G\  exp 


Ait* 

2KxKi 


dz' 


(B.7) 


&x  = 


y~v 

z 


Then 


dGx 

dz 


iAk 


2KxKiZ^ 


a*:  i  a*Gi 

+ 


del  d$i 


dal 


WJ  A 

-  -~D,(p(zf,  +  u  -  V.l,  s«,  +  0  -  V,<))iJi 


(B.8) 


If  frequency  selective  efiects  are  ignored,  then  and  By  become  parameters.  Let 


B-4 


X  -  tt 


e. 


Zt 

y-v 

Zt 


where  Zt  is  the  point  where  the  receive  or  target  is  located.  Now 


In  the  following  let  t  =0. 

The  solution  to  Equation  B.9  can  often  be  approximated  by 


HI*  t/*  tl*  V* 


gtt 

‘44  ^^“"44 


2r  -2r  -2C 

0  0  0  0  0  ^ 
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(B.IO) 


The  above  coefficients  are  easy  to  calculate.  Let  us  define  the  integral 


r/  \  \  f  zx 


(2s-2)u  zy  {zt-z) 


—  •+• 
Zt 


K)! 


(B.n) 


Then 


/(4, 0,0.0)  =  C{al) 
/(0, 4.0,0)  =  C{ol) 
/(0,0.4.0)  =  C(<7j) 
/(0.0,0.4)  C(aJ) 
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where 


=  -  In  [e"^  +  exp  (-trj)  •  (l  -  c  ^)] 

The  cross  terms  are  also  straightforward.  For  example,  let  us  calculate  C*„.  Solve 

/(a£*,0,0,O  =  C(<rJ) 

=  1  -  C(a5)/2o* 

The  other  cross  terms  are  similarly  calculated. 

Equation  B.IO  is  the  basic  equation  for  large  target,  lairge  aperture  appli¬ 
cations.  If  the  large  aperture  is  a  synthetic  aperture  radar  (SAR),  it  is  convenient  to 
place  the  u  axis  along  the  satellite  track  with  v  =  0.  Then 

G,(x,u,v)  =  G.(0,0,0)  exp  |-  1^^  +  ^  +  ^  -  -  20.,^ 

I  (B.12) 

For  generating  of  signal  structure  we  need  the  power  spectrum,  that  is,  the  Fourier 
transform  of  Equation  B.12.  Assuming  m  =  2,  let 

The  power  spectrum  is 


G,(K.,K„/f.)  = 


- No - 


-  Cl)  +  KXll  ~  Cl)  +  2KotoK.e.(C.o  +  C^C^) 


+2Ku£u‘f^ll^(^yu  +  +  C,uC^)]  | 

(B.13) 
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In  general,  the  most  stressing  case  for  a  SAR  occurs  when  the  satellite  radar  beam 
travels  perpendicular  with  respect  to  the  magnetic  field  in  the  scattering  medium.  In 
this  case 


Czy  «  Cyu  «  0 


Then 


‘ "  v/i-cL  1 - - + 1 

(B.14) 

The  y  direction  is  decoupled  from  the  x,  u  direction.  Cxu  is  less  than  one  and  relates 
to  the  thickness  of  the  scattering  medium.  In  the  limit  of  a  delta  layer,  C*u  approaches 
one  and  ifu  approaches  if*. 
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APPENDIX  C 


CALCULATION  OF  THE  “LOCAL”  PHASE  STRUCTURE 

FUNCTION 


The  “local”  structure  function  is  calculated  with  a  Fortran  call  to  CNNXM. 
The  formal  equivalence  is 

D{p{xo,yo))  =  CNNXM{/)(a;„y„),/R) 

IR  anticipates  the  need  to  iterate  over  points  along  the  propagation  line  of  sight  (LOS) 
several  times  in  calculating  some  propagation  parameters.  IABS(IR]  is  used  as  an 
index  over  LOS  integration  points.  Arrays  internal  to  CNNXM  must  be  large  enough 
to  accommodate  the  largest  number  of  integration  points.  See  conunents  inside  of 
the  CNNXM  source  code.  On  the  first  call  to  CNNXM  at  an  integration  point,  four 
variables  must  be  provided  to  CNNXM  via  the  labeled  common  block,  CNDATA.  The 
block  definition  and  the  variables  are 

REAL  N,NP 

COMMON  /CNDATA/N,NP,R,RP 
N  =  n 
NP  =  n' 


R  =  freezing  to  outer  scale  size  ratio 
RP  =  inner  to  outer  scale  size  ratio 

On  the  subsequent  calls  at  that  point,  use  the  positive  index.  The  data  in  the  labeled 
common  block  is  ignored.  Stored  data  from  the  first  call  is  used  to  provide  a  fast 
result. 


The  following  block  of  code  contains  a  basic  iteration  loop  to  calculate  a 
decorrelation  distance,  p(xo,yo)-  This  value  is  at  the  EXP(-l)  point  of  the  complex 
signal  correlation  function.  This  code  is  easily  generalized  to  multiple  integration 
points.  For  most  cases  it  is  most  efficient  to  begin  the  iteration  with  a  small  value  of 
XS=  p(xo,yo)  =RP.  The  error  in  the  final  value  for  p(xo>yo)  is  less  than  five  percent. 
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PROGRAM  TEST 
REAL  N.NP.M.LO.LF 
COMMON  /CNDATA/N.NP.R.RP 
C 

C  LO  =  OUTER  SCALE 
C 

L0=l.OE7 

C 

C  LF  -  FREEZING  SCALE 
C 

LF=1.0E5 

C 

C  P2  =  MEAN  SQUARE  PHASE  FLUCTUATION 
C 

P2=1000.0 

C 

C  SET  UP  LABELED  COMMON  BLOCK  FOR  NEGATIVE  INTEGRATION  INDEX 
C  2*N-2  =  INTERMEDIATE  SCALE  SPECTRAL  INDEX 

C  2+NP-2  *  TRANSITION  SCALE  SPECTRAL  INDEX 

C  R  =  FREEZING  TO  OUTER  SCALE  RATIO 

C  RP  =  INNER  TO  OUTER  SCALE  RATIO 

C 

N=1.6 
NP=3.0 
R=LF/LO 
RP=R/ 100.0 
C 

C  START  ITERATION 
C 

CLIM=-ALOG (EXP ( - 1 . 0) + ( 1 . 0-EXP ( - 1 . 0) ) *EXP ( -P2) ) 

XS=RP 

M=2.0 

AN0=CNNXM(XS.-1) 

XO=RP/SQRT(P2*ANO/CLIM) 

56  ANSS=CNNXM(X0.1) 

DP=P2*ANSS/CLIM 

IF(ABS(1.0-DP).LT.0.001)GO  TO  57 
M=ALOG ( ANSS/ANO) /ALOG (XO/XS) 

ANO^ANSS 
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xs*xo 

XO=XO/DP**( 1.0/M) 

IF (XO . LT . (XS/ 10.0)) X0=XS/10 . 0 
GO  TO  56 
57  CONTINUE 

C 

C  AT  THIS  POINT  THE  RESULTS  ARE 
C  M=EFFECTIVE  SPECTRAL  INDEX 

C  XO=FINAL  SOLUTION  FOR  RH0(X_0,Y_O) 

C 


The  following  code  is  the  CNNXM  Fortran  function.  The  internal  arrays 
must  be  sized  to  the  largest  number  of  LOS  integration  points  to  be  used.  The  size  is 
specified  by  the  internal  parameter  NIP. 


FUNCTION  CNNXM (X.IFLG) 

C 

C  FUNCTION  CNNXM.  VERSION  3.0,  4  APR  89 

C 

C  THIS  IS  THE  PHASE  STRUCTURE  FUNCTION  FOR  THE  WITTWER-KILB  POWER. 

C  SPECTRUM.  WHEN  USED  TO  CALCULATE  THE  SIGNAL  DECORRELATION 
C  DISTANCE,  THE  MAXIMUM  ERROR  IN  THE  DISTANCE  IS  FIVE  PERCENT 
C  THE  MAXIMUM  ERROR  IN  THE  STRUCTURE  FUNCTION,  ITSELF,  IS  ABOUT  TEN 
C  PERCENT. 

C 

C  THE  INPUT  VARIABLES  ARE: 

C 

C  LABELED  COMMON  BLOCK  /CNDATA/ 

C  (THESE  VARIABLES  ARE  USED  ONLY  WHEN  IFLG  <  0) 

C  2*N-2  =  INTERMEDIATE  SCALE  SPECTRAL  INDEX 
C  (N.GE.1.6.AND.N.LE.2.0) 

C  2*NP-2  =  SMALL  SCALE  SPECTRAL  INDEX 
C  (NP.GE.2.0.AND  NP.LE.4.0) 

C  R  =  RATIO  OF  FREEZING  TO  OUTER  SCALE 
C  (R.LE.l.AND  R.GE.l.OE-07) 

C  RP  =  RATIO  INNER  TO  OUTER  SCALE 
C  (RP.LE.R) 

C 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


FORMAL  ARGUMENTS 

X  =  DISPLACEMENT  DIVIDED  BY  OUTER  SCALE 
lABS(IFLG)  »  INDEX  OF  LOS  INTEGRATION  POINT 

( lABS ( IFLG) . LE . NIP . AND . IFLG . NE . 0) 

IFLG  <  0.  CALCULATE  AND  STORE  INTERMEDIATE  RESULTS.  THIS 
IS  USED  WHENEVER  N.  HP.  R.  OR  RP  CHANGES  FOR 
SOME  VALUE  OF  lABS(IFLG). 

IFLG  >  0.  USE  STORED  INTERMEDIATE  DATA.  ASSUMES  THAT 
N.  NP.  R.  AND  RP  RAVE  NOT  CHANGED  FOR  CURRENT 
VALUE  OF  lABS(IFLG)  SINCE  LAST  CALL. 


SAVE 

REAL  N.NP.L2 
C 

C  NIP=NUMBER  OF  LOS  INTEGRATION  POINTS 
C 

PARAMETER  (NIP*U  .E2=-7 . 1509297E-01) 

PARAMETER  (D3=6. 1788469E-01 ,C3=9 .6310179E-02) 
PARAMETER  (D4=l .3862943, C1=-6.9314718E-01) 
PARAMETER  (E7=2 . 8490703E-01 .C4»6 .6561593E-01) 
DIMENSION  XPSD(12*NIP) .FP(NIP) .E4(NIP) . JP(NIP) 
DIMENSION  E6(NIP) .L2(NIP) .IA7(NIP) . JM(NIP) 
DIMENSION  XPSDP(12*NIP) .XPSDM(12*NIP) .XN(12*NIP) 
COMMON  /CNDATA/N . NP . R . RP 
I=IABS(IFLG) 

JSTART=(I-1)*12*1 
IF ( IFLG. LT.O) THEN 
C 

C  DO  A3=d0**2 
C 


A3=R**1.2 

IF(ABS(N-NP)  .LT.O.OODTHEN 
A6=EXP(-6.4/(N*(N+6.4))) 

ELSE 

A6=((J..O+6.4/NP)/(1.0*6.4/N))**(1.0/(NP-N)) 
END  IF 

A3=A3*(1.0-A3)*A0 

C 

C  INITIAL  NORMALIZATION  CONSTANT 
C 
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A6-1 , 00+ ( (0 . 06*HP-0 . 64 ) *NP-K) . 38) *ALOG(R) 
A6»(A6*R)**(2.0*N-2.0) 

A6*l . 0/ ( 1 . 0- (HP*N) * A6/ (NP- 1 . 0) ) 

FP ( I) =0 . 24*A6* (M- 1 . 0) * ( 1 . 0+6 . 4/HP) / A3* ♦ (NP-N) 
C 

C  CONSTAHTS 
C 

L2(I)=R*R/A3 

C  X0«AL06(1 .0/SqRT(L2) ) 

X0=-0.6*AL0G*(L2(I)) 

IF(X0.LT.2.05)THEN 

A3=0.961220*XO-1.0 

XH(JSTART+l)=A3-2.3 

XN(JSTART+2)=A3-1.5 

XN(JSTART+3)=A3-1.0 

XN(JSTART+4)=A3-0.5 

XN(JSTART*6)=A3 

XN(JSTART+6)=A3+0.5 

XN(JSTART*7)»A3+1.0 

XN{JSTART*8)=A3*1.66 

JEHD'JSTARt^Q 

XN(JEN0)»A3*2.46 

ELSE 

XK(JSTART*1)=-1.36 
XN(JSTART*2)=-0.66 
XN(JSTART*3)=-0.06 
XH(JSTART*4)=0.46 
J=JSTART*6 
IF(X0.GT.2.86)THEH 
XH(J)»1.35 
XN(J*1)=X0- 1.49999 
XH(J*2)=XO*0.60 
JEND=J*5 

ELSE 

XH(J*l)=X0-0.60 

XH(J)=0.6*(XH(J*i)*XH(J-l)) 

JEHD=J*4 
END  IF 

XM(JEND-2)®X0-0.10 

XH(JEHD-1)=X0*0.46 
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XN(JEND)=X0+1.36 
END  IF 
A4=-AL0G(RP) 

807  IF(XN(JEND).LT.A4)G0  TO  812 

JEND=JEND-1 
GO  TO  807 

812  JEND=JEND+1 

XN(JEND)=A4 
I\7(I)»JEND 
XPSDM(JEND)>0.0 
JM(I)«JEND 
E4(I)=H-0.6 
E5(I)=NP-N 
XPSDP(JSTART)=0.0 
JP(I)=JSTART 
XPSD(JSTART)*0.0 
DO  900  J=JS.JEND 

A3=EXP(2.0*XN(J)) 

900  XPSD ( J) «E4 ( I ) ♦ ALOG ( 1 . 0+A3) <  E6 ( I) * ALOG ( 1 . 0^ 

1  L2(I)*A3)-XH(J) 

C 

C  2*SIH(K*X/2)**2  =  C1*(K*X)**2  .  K<=DO 

C  =  C2*K**E2  .  D0<K<=D1 

C  =  C3  .  DKK 

C 

C  Cl  *  COEFFICIENT  FOR  SMALL  (K*X)**2  =0,6 
C  D3  =  PEAK  VALUE  OF  APPROXIMATE  FUNCTION  =  1.856 
C  DO  =  K  AT  APPROXIMATE  FUNCTION  PEAK  =  SqRT(D3/(Cl*X**2)) 

C  D1  »  K  AT  END  OF  OVERSHOOT  =  D4/X.  D4  =  4.0 
C  C2  »  COEFFICIENT  FOR  INTERMEDIATE  K  =  D3/D0**E2 
C  C3  =  APPROXIMATE  FUNCTION  AT  LARGE  K  =  1,10 
C  E2  -  EXPONENT  FOR  INTERMEDIATE  K  =  AL0G(C2/03)/AL0G(Dl/DO) 
C 

C  D3=AL0G(i.866) 

C  C3»AL0G(1.10) 

C  D4»AL0G(4.O) 

C  E2»(03-D3)/(D1-D0) 

C  E2=(C3-D3)/(D4-0.6*(D3-C1)) 

C  C1=«AL0G(0.6) 

C  C4=0.6*(D3-C1) 


C  E7=E2+1.0 

C 

END  IF 
CNNXM-0.0 
IF (X.EQ. 0.0) RETURN 
XX-ALOG(ABS(X}) 

X2«X*X 

C  D0*0.6*(D3-XX-XX-C1) 

C  D1=D4-XX 

C  E2»(C3-D3)/(D1-D0) 

C  E7-1.0+E2 

D0*C4-XX 

XN ( JSTART) =AMIN1 (DO .XN ( JSTART^l) ) -6.6 
XPSD ( JSTART) =-XN ( JSTART) 

J-IA7(I) 

IF(XN(J).GT.DO)THEN 

D1=D4-XX 

C2»D3-E2*D0 

IF(XN(J).GT.D1)THEN 

JE»J-1 

JT*JE» JSTART 
DO  800  JJ* JSTART. JE 
J*JT-JJ 

800  IFCXN(J).LT.D1)G0  TO  860 

860  J=J>1 

IF(J.LT.JM(I))THEN 

JE=J.M(I)-1 

JT^JE+J 

A4=C3-XPS0(JE^l) 

F0=EXP(A4) 

DO  870  JJ^J.JE 
JS^JT-JJ 

A6«A4 

A4»C3-XPSD(JS) 

FO  *  EXP(A4) 

AO  »  A6  -  A4 

IF(ABS(AO).LT.0.0001)FK-(FN-F0)/A0 
870  XPSDM(  JS)  -XPSDMC  JS'*  1)  ♦FN*  (XN  (  JS^  1)  -XN  <  JS)  ) 

JM(I)»J 
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END  IF 

CNHXM>XPSDM(J) 

CT=EXP(D1+D1) 

CT=E4 (I) *ALOG ( 1 . 0+CT) +£6 (I) ♦ALOG( 1 . 0+L2 ( I) ) *CT) 

A6=C3+D1-CT 

F0=EXP(A6) 

A4=C3-XPSD(J) 

A1«EXP(A4) 

A0*A4-A6 

IF(ABS(AO) .LT.0.0001)A1*(A1-F0)/A0 
CNHXM«CHNXN^ A 1 * (XN ( J) >D1 ) 

ELSE 

D1=XN(J) 

CT=EXP(D1+D1) 

CT»E4(I)vAL0G(l.O^CT)^E6(I)*AL0G(l.O-»L2(I))*CT) 

A6=C2+E7*D1-CT 

F0*EXP(A6) 

END  IF 
876  J=J-1 

IF(XN(J).LT.DO)GO  TO  876 
FN=FO 

A4=C2+E2*  XN { J) -XPSD ( J) 

FO=EXP(A4) 

A0sA6-A4 

IFCABSCAO) .LT.0.0001)FN*(FN-F0)/A0 
CHNXM»CNNXM*FN* (01 -XN ( J) ) 

01*XN(J) 

A6>A4 
GO  TO  876 

876  CT»EXP(DO^DO) 

CT*E4 ( I ) *  ALOG ( 1 . O^CT) ♦£& ( I) • ALOG ( 1 . 0^L2 ( I ) ‘CT) 

A4*C2*E7»DO-CT 

A0»A6-A4 

IF(ABS ( AO) . GT . 0 . 0001 ) F0= (F0-EXP(A4) ) /AO 
CSNXM-CHNXM^FO* (Dl-DO) 

A6«Cl»3.0*D0-CT 

A4«CU2.0*XN(J)-XP$0(J) 

A1=EXP(A6> 

A0«A6-A4 

IF(ABS(AO) .GT . 0 . 0001 ) A l« ( A1 -EXP (A4) ) /AO 


C-% 


970 


CNNXM=CNNXM+X2*A1* (D0-XH( I) ) 

END  IF 

IF(J.GT.JP(I))  THEN 
JS=JP(I) 

A6«J1+2.0*XH(JS)-XPSD(JS) 

FH=EXP(A6) 

JS-JS+1 
DO  970 
A4»A6 
FO=FN 

A6=Cl+2 . 0*XN ( J J) -XPSD ( J J) 

FN*EXP(A6) 

A0=A6-A4 

IF(ABS(AO).GT. 0.0001)  FO  =  (FN-FO)/AO 
XPSDP(JJ)*XPSDP(JJ"1)+F0*(XN(JJ)-XN(JJ“1)) 
JP(I)=J 
END  IF 

CNNXM»CMNXM+X2*XPSDP ( J) 

CNNXM»FP(I)*CNNXM 

RETURN 

END 
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APPENDIX  D 


CALCULATION  OF  THE  FREQUENCY  SELECTIVE  BANDWIDTH 


This  appendix  establishes  the  basts  and  the  approximations  for  the  algo¬ 
rithm  used  to  calculate  the  frequency  selective  bandwidth,  /q.  The  transmitted  signal 
is  written  as 


stW  =  /  5(/)exp{-i2x/t}<(f  (D.l) 

y-00 

where  S{f)  is  the  Fourier  transform  of  the  waveform.  Since  is  real, 

S(f]  =  «•(-/)  . 

Abo,  the  signal  b  assumed  narrowband  so 


s(/)  =  i[P(/ +  /.)  +  >(/-/.)] 

P{f)  =  /  p(0exp{-*2x/«}  dt 

/-go 

where  /.  is  she  carrier  hettueacy  and  p(t]  is  the  basebaad  modvlation  waveform.  Also 

r  s-{mf)<if=t  . 

/-go 

The  recdved  waveform  b  written  as 

-»k(0  -  j  df  S{f)U{K)  exp  1*71  j  y/idz-  i2r/t|  (0.2) 

where  U{K)  was  defined  in  Appendix  B.  Non-essential  dependences  have  been  sup¬ 
press^  for  convenience.  Also 
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K  = 


2‘Kf 


yfe  = 


1  /|?(^)  _ .1  2jrg»ne(g) 

^  P  ^  mc»fC> 


n«(a)  =  mean  electron  density 

:  tq  =  classical  electron  radius  (2.82  x  10 


me* 


/p  =  plasma  frequency 


After  a  bit  of  algebraic  manipulation,  the  signal  time  of  arrival  moments  can  be  written 
as 


("4  (<)<«  =  <Vt  “P  /  \/«j  j 

+  exp  J  y/eidz^  ( j 


(G-(KuK,}S{fuap 


(D.3) 


where 


GiK^,Ki)=^U-iKi)UiKi) 

Since  the  signal  is  narrowband,  the  integrand  in  Equation  D.3  is  nonzero 
and  for  values  of  /j  and  fi  near  i/c  and  Kj  and  Ki  near  Thus  we  will  develop 

a  solution  for  r(iritK:)  as  an  expansion  in  terms  of 
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Since 


_  Kt-Kt  h-h 


<PG  /  iuV'i'a 
df"^  \  Kic)  dif* 

we  will  be  able  to  evaluate  terms  in  Equation  D.3.  The  final  result  required  is  the 
contribution  to  the  standard  deviation  of  the  energy  time  of  arrival  from  the  angular 
scattering  of  the  signal  energy.  Let. 

r  =  /**  di 

J-OO 


The  desired  result  is 


~  ~  ^  iCfctUr  (^•^) 

It  should  be  noted  that  Equation  D.3  is  a  very  general  result.  Various  terms  account 
for  delay  due  to  the  angular  scatter,  the  finite  length  of  the  waveform,  and  the  mean 
total  electron  content  effects. 

The  equation  for  the  G{KuK^)  is 


dz 


d*G\  zdG 
^  j  z  dx 


z  dy 


i£ik  fd*G 
KiKi  Vdz» 

1  [if?  +  Kll  dal.  ^  ,  ...  ^ 

-2  1-^1  +  df  I*  -  ® 


where  Ak  =  Ki~  Ki  and  D^(p(z,y))  is  the  phase  structure  function. 


Let 


(D.5) 


I.  = 


z 

z 


6, 


y 

z 
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Equation  D.5  can  now  be  written 


dz  ^\d$i^  del)  \ 


i  +  n  +  7ir  + 


dal 


+  (l  +  q  +  q*  +  •  •  •)  (1  ~  ■D^(p(2^x»  ^ 


dal 

where  is  evaluated  at  Jf  =  Jfj. 


Let 


Then 


Now  let 


Then 


i[M) 

1  —  D^{p{z$xi  ^^v)) 


=  [  dzf{z) 

Jo 

=  1  +  Axel  +  A^el  +  Axy^x^y 

+  +  By^J  +  Bsy6*dJ  4-  .  .  . 

=  »'exp|-/[(l  +  ,  +  ?.,’)^]| 


dz 


j  ,  fd*W 

del ) 

da^ 

+(!  +  »?  +  »)’  +  ••  (l  +  AxBl  +  •  *  •) 


iV=:IVo  +  qiV,+?^*Wa 


(D.6) 


(D.7) 
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H'o  =  «apM^(l  +  ^«^  +  •■■)|}  (D.8) 

Also 

I  -  ^  ■  = +§  (i + P-9) 

Assume 

Wi^WoTsMi^' 

JQ 

then  it  is  clear  that 

j  2 

91  Wivo  =  (i + ^*5 + •  •  •)  ‘*'1  P-io) 

Also 


Vi  Wo  = 


■)  w'o  +  7’  ^|;(i +  '»■«!  + 

(iff* 

+/  (l  +  + 


Finally, 


(D,ll) 

i  idoi 
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since  in  the  end  we  let  6z  =  ^y  —  0. 


ForWj 


Wt  =  Wog,{z) 


2K,z^ 


dffl 


+  -^  (l  +  AJl  +  •'■)  {Wo  +  Wi)  (D.13) 


where 


W2=Wo  rgi{z')dz’ 

Jo 


After  much  algebra 


gt{z)  = 


2K2Z^ 


^(2^.+2X,)]-^ 


§  (2A.  +  2.1,)]  {2A.  +  2A.))  +  ^i]  + 

^  [-2^ 

+/  (24B.  +  +  24B»)  j  j  +  ^  (2A.  +  2A„)j  | 


Finally 


(D.14) 


G-exp  |-/ +  ^1?*  +  •  •  •)  1  [1  +  fll(gi}  +  +  •  •  •] 
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After  expanding  the  exponential 


G  =  1  +  rjGi  +  ri^Gi  + 


(D.15) 


(D.16) 


Gi{z)  = 


•'))l  I 


+2J' 


‘  ["2^^  (■^ 


(24J3js  +  ^Bxfi  +  24B| 


dal 

dz 


(D.17) 


Moments  in  Equation  D.3  can  now  be  evaluated.  Let  us  first  calculate  the  contribution 
to  the  first  moment  from  the  mean  plasma  density.  This  requires  the  calculation  of 

^  (exp  {-.'Jf  £  Vide})  =  ^  £  s:(ij<k|  exp  {-iK £  Vidz] 


Thus 
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+  '•<(»)<<»)  (D.lSa) 


=  7  + 

The  first  term  is  the  signal  propagation  time  at  the  speed  of  light  and  the  second  term 
is  the  well  known  group  delay. 

Now  let  us  look  at  the  contribution  to  the  mean  delay  from  the  scintillation 

term 


(D.19a) 

(D.19b) 

(D.20) 

After  a  horrible  amount  of  algebra  (the  proverbial  exercise  for  the  reader),  the  variance 
of  the  delay  from  the  scatter  is 

ol\a 

(2x/,)*  2(27r/,)*'  1  [2* 

• 

/  (24B.  +  +  24B,)  j  +  2/*  (^2^) 

(D.21) 

The  first  term  is  from  the  group  delay  jitter  and  we  drop  it  since  creates  no  v/aveform 
distortion.  Group  delay  jitter  is  generally  handled  separately.  Now  let  us  look  at  the 
B  terms.  For  the  often  used  “if-*”  spectrum 


1 

2t^Ll 


24B, 

A. 

The  ratio  of  the  contribution  from  these  terms  is  approximately 

jt 


ln(L«/x) 

2L\ 


where 


L»  =  outer  scale 
4  =  inner  scale 

This  term  is  always  small  for  satellite  frequencies  and  applications.  Thus  it  is  dropped 
for  efficiency.  The  B  term  also  lead  to  pecular  results.  As  4  goes  to  zero,  the  delay 
becomes  infinite.  The  signal  energy,  however,  is  still  reasonably  described  by  the  A 
terms.  The  divergence  is  traceable  to  a  common  problem  of  moment  methods  with 
functions  that  do  not  go  to  zero  rapidly  for  large  arguments.  Sufficiently  high  moments 
will  diverge  while  the  function  is  still  integrable.  Thus  dropping  the  B  terms  provides 
insurance  against  unreasonable  results. 

For  the  power  spectrums  described  earlier  in  this  report  the  A  coefficients 
are  related  to  the  electron  density  scale  sizes. 


A. 

A„ 


-ft 


L*  L’  -  L* 

"tv 


Llz^ 


LlLl  -  LI 


*v 


-B, 


2LxyZ* 

LlLl  -  L%, 
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The  derivation  of  B%  is  described  earlier  leading  up  to  Equation  1.21.  This  particular 
method  insures  reasonable  results  even  if  the  second  moment  in  delay  is  not  well 
defined.  Finally  we  set 


f'*  = 

'  4ir«<r>|o 

k  =  + 


where  0  <  Ci  <  0.25.  The  Ci  parameter  was  added  to  ease  the  numerical  handling 
of  the  generalized  power  spectrum  in  the  CIRF  computer  program  which  provided 
explicit  representations  of  scintillated  signals.  Currently  it  is  chosen  in  ACIRF,  the 
successor  to  CIRF,  depending  on  the  specific  channel  variation  being  run. 

The  integral  in  Equation  D.21  can  be  simplified.  The  original  form  can  be 
expressed  as 


*•  ^ 


dz  z^f(z)  —  dx{zt  —  a;)*/(x)  j  dz  z^f{z) 


The  latter  form  is  much  easier  to  implement  and  is  strongly  recommended. 
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APPENDIX  E 


AN  ALGORITHM  FOR  CALCULATING  THE  SCINTILLATION  AND 
ABSORPTION  OF  RADIO  SIGNALS 


This  appendix  lists  SUBROUTINE  PROP  which  calculates  the  scintilla¬ 
tion  and  absorption  of  radio  signals.  The  inputs  and  outputs  are  documented  in  the 
comments  isection  at  the  beginning  of  the  subroutine.  This  routine  implements  the 
algorithm  described  in  the  body  of  this  paper  for  an  environment  consisting  of  a  series 
of  slabs  each  with  its  particular  properties.  Thus  the  subroutine  can  handle  most 
situations. 


SUBROUTINE  PROP(MOOE) 

C 

C  PROP.  VERSION  6.  10  FEB  90 
C 

C  THIS  ROUTINE  CALCULATES  PARAMETERS  DESCRIBING  THE  PROPERTIES 
C  OF  RADIO  SIGNALS  THAT  RAVE  PROPAGATED  THROUGH  AH  ABSORBING 
C  AND  STRUCTURED  MEDIUM.  THE  ABSORPTION  MODEL  INCLUDES 
C  ELECTRON- ION  AND  ELECTRON-NEUTRAL  COLLISIONS.  THE  PARAMETERS 
C  DESCRIBING  THE  SCINTILUTION  OF  RADIO  SIGNALS  ASSUMES 
C  PROPAGATION  THROUGH  A  STRUCTURED  PUSMA  DESCRIBED  BY  THE 
C  WITTNER/KILB  STRUCTURE  MODEL.  THE  LIMITS  ON  INPUT  PARAMETERS 
C  ARE  DESCRIBED  BELOW.  FOR  EFFICIENCY.  NO  CHECKS  ARE  MADE  TO 
C  TEST  THE  INPUTS  FOR  LEGALITY  EXCEPT  FOR  THE  HUMBER  OF 
C  ENVIRONMENT  UYERS. 

C 

C  THIS  PROGRAM  ASSUMES  THAT  ANY  INPUTS  HAVE  AT  LEAST  A  REMOTE 
C  RESEMBLENCE  TO  PHYSICAL  REALITY.  SINCE  I  CAN  NOT  PROTECT 
C  AGAINST  ALL  CONCEIVABLE  POSSIBILITIES.  I  ONLY  CLAIM  THAT  THE 
C  REASONABLENESS  OF  THE  OUTPUTS  IS  AT  LEAST  COMPARABLE  TO  THE 
C  REASONABLENESS  OF  THE  INPUTS. 

C 

C  THIS  SUBROUTINE  WAS  PROGRAMMED  BY  DR.  LEON  A.  WITTWER. 

C  QUESTIONS  HAY  BE  REFERRED  TO  DR.  WITTWER  AT  703-326-7028. 

C 

C  EBRORS  OR  WARNINGS  ARE  OUTPUT  VIA  WRITE(« . )  CALLS. 
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C  IS  ASSUMED  DEFAULTED  TO  THE  TERMINAL. 

C 

C  MODE  CONTROLS  THE  PROPERTIES  CALCUUTED  AND  RETURNED.  SEE 
C  BELOW  FOR  DETAILS. 

C 

C  LABELED  COtOION  /PRINPT/  CONTAINS  INPUT  TO  THE  SUBROUTINE 
C 

C  F  «  CARRIER  FREQUENCY (HZ) . 

C  NLYRS  >  NUMBER  OF  ENVIRONMENT  UYERS. 

C  (NLYRS. LE. NIP  WHERE  NIP  -  MAXIMUM  NUMBER  OF  LAYERS 

C  ALONG  THE  LINE  OF  SIGHT.  NIP  IS  SET  IN  A  PARAMETER 

C  STATEMENT  BELOW.) 

C  LQSV(I)  =  VECTOR  PAR.ALLEL  TO  LIKE  OF  SIGHT.  THE  USER 
C  COORDINATE  SYSTEM  MUST  BE  A  RIGHT  HANDED 

C  CARTESIAN  SYSTEM.  THIS  VECTOR  DEFINES  THE  LOCAL 

C  Z  COORDINATE  ALONG  THE  LINE  OF  SIGHT. 

C  LOSd)  =  LINE  OF  SIGHT  COORDINATES  iKM)  OF  ENVJHCi.MENT 
C  UYER  CENTERS.  ONE  END  OF  THE  LINK  IS  AT  Z=0 

C  AND  THE  OTHER  END  IS  AT  Z-  L0S(NLYRS<^1) .  ANTENNA 

C  REUTED  OUTPUTS  (CPT,  CQT.  AND  TWPCC)  ASSIWE  THAT 

C  THE  ANTENNA  IS  AT  2=L0S(HLBYS*S) . 

C  (LOSd+D.GT.LOSd)  FOR  ALL  I) 

C  DLOSd)  «  LAYER  THICKHESS(KM)  OF  UYER  CENTERED  AT  LCSd) . 

C  NEd)  -  MEAN  ELECTRON  DENSITY (/CM**3)  IN  .  'YER  CSKTEREr 
C  AT  LOSd) . 

C  HE2(I)  s  ELECTRON  DENSITY  VARIAKa(/CM*^6)  IN  UYLR 
CENTERED  AT  LOS(I). 

NDd)  -  NEUTRAL  DEHSITY(/CM**3)  IH  LAYER  CENTERED  AT  LOSd). 
TEd)  »  MEAN  PL  ASMA  TEMPERATURE  (DEG  KELVIN)  IN  UYER  CENTERED 
AT  LOSd) . 

BV(J.I)  -  VECTOR  PARALLEL  TO  THE  EARTH'S  54AGNETIC  FIELD 
AT  LOSd),  THE  MODEL  ALLOWS  FOR  HONISOTROPIC 
STRUCTURE  IH  ONE  DlkfiCTIOH  DESIGNATED  BY  BV(J.I). 

J  IS  THE  CARTESIAN  COORDINATE  INDEX.  STRUCTURE 
IS  ASSUMED  ISOTROPIC  ABC UT  THE  VECTOR  BV. 

C  LOd)  »  OUTER  SCALE(Kk,  PERPE?  OICULAR  TO  BV  AT  LOS(I). 

C  LFd)  »  FREEZING  SCALE(KM)  PE''^ENDICULAR  TO  BV  AT  LOS(I). 

C  (LF(I).LE.L0(I).AND.LF(I).GE.l.OE-7*L0(I)) 

C  Lid)  -  INNER  SCALE(KM)  PERPENDICULAR  TO  BV  AT  LOS(I). 

C  aid).LT.LF(I)) 
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C  LT(I)  •  CUTER  SCALE(KM)  PARALLEL  TO  BV  AT  LOS (I) .  THE 
C  FREEZING  SCALE  PARALLEL  TO  BV  VICTOR  IS  LT(I)* 

LF(I)/LO(I).  THE  INNER  SCALE  PARALLEL  TO  3V  IS 
LT(I)*LI(I)/LO(I). 

(LT(i;.LE.lOO.O*L0(I)) 

N(I)  2.0*N(I)*2.0  IS  INTER34EDIATE  SCALE  SPECTRAL  INDEX. 

(N(I).GE.1.6.AND.N(I).LE.2.0) 

NP(I)  2.0«KP(I)-2.0  IS  SMALL  SCALE  SPECTRAL  INDEX. 

(NP(I).GE.2.0.AND.NP(I).LE.4.0) 

VST(J.I)  »  VECTOR  STRUCTURE  VELOCITY (KM/SEC)  AT  LOS(I). 

J  IS  THE  CARTESIAN  COORDINATE  INDEX. 

VTR(J)  -  TRANSMITTER  VELOCITY (KM/ SEC) .  J  IS  THE  CARTESIAN 
COORDINATE  INDEX. 

VRE(J)  *  RECEIVER  VELOCITY (KM/SEC) .  J  IS  THE  CARTESIAN 
COORDINATE  INDEX. 

THE  OUTPUT  IS  DETERMINED  BY  THE  VALUE  OF  MODE. 


C 

C 

C 

G 

C 

C 

C 

C 

C 

C 

C 

C 


MODE  >  0 

LABELED  COMMON  /PROUTO/ 

TO  »  SIGNAL  DECORRELATION  TIMECSEC) 

THPCC  •  TWO  WAY  PATH  SIGNAL  CORRELATION  COEFFICIENT 
ASSUMING  THAT  THE  TRANSMITTSR/RECEIVER  IS  AT 
Z-LOS(NLYRS^l) . 

CPT  -  TIME  SPACE  GROSS  CORSEUTION  COEFFICIENT  IN  P 
DIRECTION  AT  Z*LOS(HLYRS*l) 

GOT  -  TIME  SPACE  CROSS  CORRELATION  COEFFICIENT  IN  Q 
DIRECTION  AT  Z*LOS(HLYRS*l) 

FO  >  FREGUENCY  SELECTIVE  BANDWIDTH (HZ) 

X2  MEAN  SQUARE  LOG  AMPLITUDE  FLUCTUATION 
KA  >  ABSORPTION (DB) 

AT  -  ANTENNA  TEMPERATURE(DEG  KELVIN)  AT  RLJEIVER. 

ASSUMES  FIREBALL  FILLS  ENTIRE  ANTENNA.  THE 
RETURNED  TEMPERATURE  IS  FOR  THE  PLASMA  BETWEEN 
THE  TRANSIITTER  AND  RECEIVER.  IT  DOES  NOT 
INCLUDE  ANY  PLASMA  BEHIND  THE  TRAN^ITTER. 

LPV  =  VECTOR  IN  DIRECTION  OF  MINIMUM  DECORRELATION 
DISTANCE  AT  Z^LOSCNLYRS*!)  DEFINED  AS  THE  P 
DIRECTION. 

LP  »  MINIMUM  DECORREUTION  DISTANCE  AT  Z»LOS(NLYRS^l)(M) . 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


LQ  >  MAXIMUM  DECORRELATIOH  })ISTAHCE  AT  Z>=LOS(NLYRS^l)  (M) . 
THE  DIRECTION  IS  PERPENDICULAR  TO  THE  VECTOR  LPV 
AND  THE  LOS. 

LPPV  -  VECTOR  IN  DIRECTION  OF  MINIMUM  DECORRELATION 
DISTANCE  AT  Z»0. 

LPP  =  MINIMUM  DECORRELATION  DISTANCE  AT  Z=0(M) . 

LQP  ==  MAXIMUM  DECORRELATION  DISTANCE  AT  Z=0(M) . 

THE  DIRECTION  IS  PERPENDICULAR  TO  THE  VECTOR  LPPV 
AND  THE  LOS. 

CPPP  =  TRANSMITTER-RECEIVER  CORRELATION  COEFFICIENT  IN 
DIRECTION  OF  MINIMtJM  DECORRELATION  DISTANCE  AT 
Z=0. 

MODE  =  1 

ALL  OF  THE  ABOVE  OUTPUTS  AND 
LABELED  COMMON  /PROUTl/ 

P2  »  MEAN  SQUARE  PHASE  FLUCTUATION (RAD* *2) 

TEC  =  TOTAL  ELECTRON  C0NTENT(/CM**2) 

LOEF  »  EFFECTIVE  OUTER  SCALE (KM) 

LFEF  =  EFFECTIVE  FREEZING  SCALE (KM) 

LIEF  »  EFFECTIVE  INNER  SCALE(KM) 

NEF  2.0*NEF-2.0  IS  EFFECTIVE  INTERMEDIATE  SCALE 
SPECTRAL  INDEX 

NPEF  2.0*NPEF-2.0  IS  EFFECTIVE  SMALL  SCALE  SPECTRAL 
INDEX 

VEF  =  EFFECTIVE  TEC  VELOCITY 

17  VEF  IS  ZERO  THEN  LOEF.  LFEF.  LIEF,  NEF,  NPEF. 

KR.  IR.  DLM(J).  AND  DPM(J)  ARE  UNDEFINED. 

MODE  =  2 

ALL  OF  THE  ABOVE  OUTPUTS  AND 
LABELED  COMMON  /PR0UT2/ 

KR  =  RAYLEIGH  WAVENUMBERd  .0/KI4) 

FR  =  RAYLEIGH  FREQUENCY (HZ) 

DL  =  MEAN  GROUP  DELAY (SEC) 

DLM(l)  =  3  SIGMA  DELAY(SEC) 

DLM(2)  =  3  SIGMA  DELAY  RATE(SEC/SEC) 

DLM(3)  *=  3  SIGMA  DELAY  ACCELERATI0H(SEC/SEC**2) 

DLM(4)  =  3  SIGMA  DELAY  JETJC (SEC/SEC* *3) 

DP  =  MEAN  PHASE (R\D) 

DPM(l)  =  3  SIGMA  PHASE(RAD) 

DPM(2)  =  3  SIGMA  DOPPLER(HZ) 
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C  DPM(3)  »  3  SIGMA  DOPPLER  RATE(HZ/SES) 

C  DPM(4)  «  3  SIGMA  JERK (HZ/SEC* *2) 

C 

C  SET  PARAMETERS 
C 

C  NIP  >  MAXIMUM  HUMBER  OF  ENVIRONMENT  LAYERS.  NIP  MUST  ALSO 
C  BE  SET  IN  FUNCTIONS  CNNXMI  AND  CNNXM. 

C  RC  -  RAYLEIGH  CRITERIA 

C  PND  «  DEFAULT  PHASE  NOISE 

C  Cl  -  FREQUENCY  SELECTIVE  PARAMETER 
C  XLEST  -  INITIAL  SCALE  ESTIMATE 
C  PI02  »  PI/2 

C  PI04  =  PI/4 

C  TMPK  METERS/KILOMETER  =  1000,0 

G 

PARAMETER  {NIP“1 1 . RC=0 . 1 . PKD»0 . 025 , C1«0 . 25 , XLEST=1 . OE-5) 
PARAMETER  (TMPK=1000 .0 .PI02=1 . 5707963 .PI04=0 . 78539816) 

C 

REAL  LOS.NE.LO,LF,LI.LT.N.NP.NE2.LOSV.ND 

RFJiL  LP.LPV.LQ.KA 

REAL  LOKF.LFEF.NEF.NPEF.LIEF 

REAL  LPP.LPPV.LQP 

REAL  KR 

BEAL  LU.LV.LY2.LUP.LVP 
*»EAL  LOSVA 

DIMENSION  LOSCNIP^l)  .DLOS(NIP)  .NE(NIP)  .TE(NIP)  .BVO.NIP) 
DIMENSION  LFCNIP) .LI(NIP) .LT(HIP) .N(NIP) ,HP(NIP) ,NE2(NIP) 
DIMENSION  VSTO.NIP)  .LO(NIP)  ,L0SV(3)  .KD(NIP) .  VTR(3)  .VRE(3) 
DIMENSION  LPV(3) 

DIMENSION  LPPV(3) 

DIMENSION  DLM(4).DPM(4) 

DIMENSION  LY2(N1P) .DP2(NIP) ,VT(NIP) .UT(NIP) .UVT(HIP) 
DIMENSION  X(3) .ARGCNIP) .U(3) .V(3) .C2(NIP) .XN(20) 

DIMENSION  UV(NIP).VV(HIP) 

COMMON  /PRINPT/ F . NLYRS . LOSV . LOS . OLOS . JE . NE2 . ND . TE .BV . LO . LF 
1  .LI.LT.N.NP.VST.VTR.VRE 

COMMON  /PRQUTO/TO , TWPCC . FO . X2 . KA . AT . LPV . LP . LQ . CPT . CQT . 

1  LPP.LPPV.LQP.CPPP 

COMMON  /PROUTl /P2 . TEC . LOEF . LFEF . LIEF . NEF . NPEF . VEF 
COMMON  /PR0UT2/KR.FR.DL.DU4.0P.DPM 
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COMMON  /CNDATA/QN.QNP.R.RP 
COMMON  /CNNCOM/CLIM . ARC .DP2 
C 

C  INITIAL  SET  UP  FOR  PROPAGATION  CALCULATION  AND  CALCULATE 
SIMPLE  INTEGRAL  PROPERTIES (P2.KA. TEC. AT) . 

IF(NLYRS.GT.NIP)THEN 
1fRITE(*.6) 

6  FORMATCIX.TATAL  ERROR—TO  MANY  LOS  INTEGRATION*. 

1  *  POINTS*) 

PAUSE  *CR  TO  END* 

STOP 
END  IF 
C 

CL-  INDEX  OF  LARGEST  LOSV  COMPONENT 
C 

TPl-ABS(LOSV(l)) 

L-1 

DO  30  >2.3 

IF(ABS(L0SV(J)).LT.TP1)G0  TO  30 
L-J 

TPl»ABS(LOSV(J)) 

30  CONTINUE 

ZT-L0S(NLYRS+1) 

LOSVA-SQRT (LOSV ( 1 ) *LOSV ( I ) ♦LOSV (2) •LOSV (2) ♦LOSV ($) *LOSV (3) ) 

P2«0,0 

KA-0.0 

AT-0.0 

TEC-0.0 

VZ-0.0 

X2-0.0 

DO  1000  I-l.NLYRS 

TP3-BV(l.I)*BV(i.I)+BV<2.I)«BV(2.I)^BV{3.I)*BV(3.I) 

C 

C  CALCUUTE  LOCAL  S  DIREcTFON  VECTOR 
C 

X(l)»BV(2a‘i»LS^«VC3)  •EV<3.I)*L0Sy(2) 

X(2)*BV(3 .  X>  *L0SV<i)  Cl .  I)  •lOSVO) 

X(3)-BV(1.2;  *lQSVC2)-dV(3.I)*L0SV(l) 
TPI-X(1)*X(1>^X(2>»X(5b)^vC3)*X(3) 


IF(TP1/(TP3*L0SVA*L0SVA) .LT.l .0E-08)THEH 
LY2(I)-L0(I)*L0(I) 

J=L*1 

IF(J.0T.3)J=1 

X(J)*1.0 

X(L)=-LOSV(J)/LOSV(L) 

J=L-1 

IF(J.EQ.0)J=3 

X(J)*0.0 

TPl*1.0+X(L)*Xa) 

ELSE 

C 

C  CALCULATE  LOCAL  SQUARE  OF  OUTER  SCALE.  LY2.  IN  LOCAL  Y 
C  DIRECTION.  THE  LOCAL  SQUARE  OF  OUTER  SCALE  IN  X  DIRECTION 
C  IS  LO*LO. 

C 

TP4= (BV ( 1 . 1) ♦LOSV ( 1 ) +BV (2 . 1) *LOSV (2) ♦BV (3 . 1) *LOS  V  (3) ) 
LY2(I)*((L0(I)*TP4)**2+LT(I)*LT(I)*TP1)/(L0SVA*L0SVA 
1  *TP3) 

END  IF 

TP2»SQRT(TP1) 

IFCI.EQ.DTHEN 

C 

C  SET  INTERNAL  REFERENCE  VECTORS.  U  AND  V 
C 

U(1)=X(1)/TP2 

U(2)=X(2)/TP2 

U(3)=X(3)/TP2 

V(1)=(L0SV(2)*U(3)-L0SV(3)*U(2:>/L0SVA 
V(2)=a0SV(3)*U(l)-L0SV(l)*U(3))/L0SVA 
V(3)-(L0SV(1)*U(2)-L0SV(2)*U(1))/L0SVA 
END  IF 
C 

C  SET  UP  LOCAL  TRANSFORMATION  MATRICES  FOR  U.V  COORDINATES 
C 

J=L^l 

IF(J.GT.3)J=1 

M*J^1 

IF(M.GT.3)M=1 

TP1»  (X  ( 1)  *0  ( 1 )  ♦X  (2)  *0  (2)  ♦X  (3)  *0(3)  )  /TP2 
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TP2-L0SVA* (X< J) *U(M) -X(M) *U<J) )/(TP2*L0SV(L) ) 
TP3-L0(I)*L0(I) 

UT(I)-TP1*TP1/TP3+TP2*TP2/LY2(I) 

VT(I) =TP1*TP1/LY2 (I) ♦TP2*TP2/TP3 

UVT ( I) -2 . 0*TP 1 *TP2* ( 1 . 0/TP3-1 . 0/LY2 ( I) ) 

C 

C  CALCUUTE  LOCAL  VELOCITIES  IN  U  AND  V  DIRECTIONS 
C 

DO  901  >1.3 

901  X(J)-VST( J . I) - (LOS(I) *VRE( J)+(ZT-LOS(I) ) *VTR( J) ) /ZT 

UV(I)»U(1)*X(1)+U(2)*X(2)+U(3)*X(3) 
VV(I)»V(1)*X(1)+V(2)*X(2)+V(3)*X(3) 

C 

C  INITIALIZE  CNNXM  FOR  ITH  LAYER 
C 

C  QN  =  n 
C  QNP  =  n‘ 

C  R  »  FREEZING  SCALE  /  OUTER  SCALE 
C  RP  »  INNER  SCALE  /  OUTER  SCALE 
C 

QN=N(I) 

QNP'NPd) 

R»LF(I)/LO(I) 

RP=LI(I)/LO(I) 

TPZ-CNNXMCRP.-I) 

C 

C  TPl=N2(n.n’.R) 

C 

TPl=1.0*((0.06*qHP-0.64)*QNP>0.38)*ALOG(R) 

Q2(I)=1.0/(1.0-(QNP-QN)*(TPl*R)**(QN^QN-2.0)/(QNP-1.0)) 

C 

C  TP2=N3(n.n*.R) 

C 

TP2»(0 . 12*QNP- 1 . 3) *QNP>2 . 18^ ( (-0 , 1*QNP+ 1 . 1 1) *QHP-0 . 49) 

1  /SQRT(R) 

TP2»TP2/SQRT(1 .0>0.3/R) 

TP2*! . 0/ ( 1 .0- (QNP-QH) * (TP2  *R) (QN^QN-3 .0) / (QNP- 1 . 6) ) 

C 

C  TPS  »  EFFECTIVE  PARALLEL  SCALE  SIZE. 

C 
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TP3»L0(I)*LT(I)/SQRTaY2(I)) 


C 

C  THE  FOLLOWING  IS  USED  TO  ELIMINATE  EXPLICIT  CALCULATION  OF 
C  THE  GAMMA  FUNCTION 

C  (0 . 12+0 . 77/N) «GAMMA(N-0 . 5) / (SqRT(PI) *GAMMA(N) ) 

C 

C  TPl  »  INTEGRATED  VARIANCE  OF  PROCESS  DESCRIBED  BY  NE2 
C 

TPl-2 . 0* (QN- 1 . 5) *TP2*TP3*DL0S (I) / (Q2 (I) * (QN- 1 . 0) * (0 . 12 
1  ♦0.77/QN)) 

C 

C  DP2(I)  «  INTEGRATED  PHASE  VARIANCE  FOR  ITH  UYER 
C 

DP2(I)=7.12E+06*TP1*NE2(I)/(F*F) 

P2=P2+DP2(I) 

C 

C  GET  DOMINANT  LAYER.  TRY  TO  SAVE  SOME  TIME 
C 

TP1»DP2(I) * (UT(I) ♦VT(I) ) 

IF(TP1.GT.X2)THEN 
IDL=I 
X2*TP1 
END  IF 
C 

C  AT  »  ANTENNA  TEMPERATURE (DEG  KELVIN) 

C  KA  =  ABSORPTION (DB) 

C 

C  TP16  =  ELECTRON  ION  ABSORPTION (DB/KM) 

C 

TP1^SQRT(NE2(I)+NE(I)*NE(I)) 

TP2=1.8*TP1*AL0G(1.3E16*TE(I)**3/(F+F))/TE(I)**1.5 

TP16=^4.6E4*TPl*TP2/(39.44*F*F+TP2*TP2) 

C 

C  TP14  =  ELECTRON  NEUTRAL  ABSORPTION (DB/KM) 

C 

TP3*2.0E-11*ND(I)*TE(I) 

TP4» ( 157 . 0+F/TP3) • +0 . 78 
TP4»2 . 0+ ( 1 . 0-TP4 ) / ( 1 . 0+TP4 ) 

TPl6=TP3/(8.8+F) 

TP3«TP3* (0 . 8+0 . 2* ( 1 . 0-TP16) / ( 1 . 0+TP16) ) 
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TP14«4 . 6E4*ME (I) ♦TP3/ ( (6 . 3*F*TP4) **2+TP3*TP3) 
TP2«TP16+TP14 
TP1»EXP (-0 . 23*TP2*DL0S ( I) ) 
AT»TE(I)*(1.0-TP1)+AT*TP1 
KA*KA+TP2*DL0S(I) 

TEC  *  TOTAL  ELECTRON  CONTENT (/CM* *2) 

TEC«TEC+1.0E6*NE(I)*DL0S(I) 

ACCUMUUTE  FOR  EFFECTIVE  PUSMA  LOS  VELOCITY 

VZ»VZ+ (VST ( 1 . I) *LOSV ( 1 ) +VST (2 . I ) *LOSV (2 ) +VST (3.1) 
1  *L0SV(3))*DP2(I) 

1000  CONTINUE 
C 

C  VZ  »  EFFECTIVE  PLASMA  VELOCITY  PARALLEL  TO  LOS 
C 

IF(P2.GT.0.0)THEN 

VZ=VZ/(L0SVA*P2) 

ELSE 

VZ=0.0 
END  IF 


C 

C  BEGIN  CALCULATION  OF  PROPERTIES  THAT  REQUIRE  ITERATION  OVER 
C  LOS  INTEGRALS  OF  THE  STRUCTURE  FUNCTION  AND  OTHER  COMPLICATED 
C  THINGS 
C 

C  CALCUUTE  LU.  LV.  AND  CUV 
C 

DO  1100  I=1.NLYBS 
1100  ARG(I)*SQRT(UT(I))*LOS(I)/ZT 
LU«CHNXMI (NLYRS . XLEST) 

DO  1200  1=1, NLYRS 
1200  ARG (I) •SQRT ( VT ( I ) ) ‘LOS ( I ) /ZT 
LV-CHNXMI (HLYRS . XLEST) 

TP6»1.0 

IF(UVT(IDL) .LT.0.0)TP6=-TP6 
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X2  -  0.56 

1326  DO  1300  I=1.NLYRS 

1300  ARG(I)«SQRT(LU*LU*UT(I)+LV*LV*VT(I)-TP6*LU*LV*UVT(I))* 
1  LOS(I)/ZT 
CUy=CKNXMI(KLYRS.0.1) 

IF(CUV.GT.X2>G0  TO  1320 
TP6=-TP6 
X2  -  0.0 
GO  TO  1325 

1320  CUV=TP6*(1.0-CLIM/(2.0*CUV*CUV)) 

C 

C  CALCUUTE  LP.  LPV,  AND  LQ 
C 

TP4»LU/LV-LV/LU 
IF(ABS(TP4) .GT.1.0E-7)THEN 
TP2=ATAN (2 . 0*CUV/TP4) /2 .0 

ELSE 

TP2=PI04 

ENDIF 

TP20=SIN(TP2) 

TP21«C0S(TP2) 

TP1=TP20/LV 

TP2*TP21/LU 

TP22  =  TP1*TP1+TP2*TP2-2.0*CUV*TP1*TP2 

TP1-TP20/LU 

TP2=TP21/LV 

TP23  =  TP1*TP1+TP2*TP2+2.0*CUV*TP1*TP2 
IF(TP22  .GT.  TP23)THEN 

IF(TP23 . LT . 1 . 0E-04*TP22)TP23=1 .0E"04*TP22 
LP*1.0/SQRT(TP22) 

LQ»1.0/SqRT(TP23) 

ELSE 

IF {TP22 . LT . 1 . 0E-O4*TP23)TP22» 1 . 0E-04*TP23 
LP*1.0/SQRT(TP23) 

LQ=1.0/SQRT(TP22) 

TP1=TP20 

TP20=TP21 

TP21=-TP1 

ENDIF 

DO  1616  J^l.S 
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1615  LPV(J)=U(J)*TP21+V(J)*TP20 
C 

C  CALCULATE  LUP,  LVP,  AND  CUVP 
C 

DO  3100  I=1.NLYRS 

3100  ARG(I)=SQRT(UT(I))*(ZT-LOS(I))/ZT 
LUP=CNNXMI (NLYRS . XLEST) 

DO  3200  I=1.HLYRS 

3200  ARG(I) «SQRT(VT(I) ) * (ZT-LOS (I) ) /ZT 
LVP»CNNXMI (NLYRS , XLEST) 

X2»0.B6 

3325  DO  3300  1=1, NLYRS 

3300  ARG(I)=SqRT(LUP*LUP*UT(I)+LVP*LVP*VT(I) -TP6*LUP*LVP*UVT(I) 
1  )*(ZT-LOS(I))/ZT 
CUVP=CNNXMI (NLYRS . 0 . 1 ) 

IF(CUVP.GT.X2)G0  TO  3320 

TP6=“TP6 

X2=0.0 

GO  TO  3325 

3320  CUVP=TP6* (1 .0-CLIM/ (2 .0*CUVP*CUVP) ) 

C 

C  CALCULATE  LPP.  LPPV,  AND  LQP 
C 

TP4=LUP/LVP-LVP/LUP 
IF(ABS(TP4) .GT.1.0E-7)THEN 
TP2=ATAN(2 .0*CUVP/TP4) /2 .0 

ELSE 

TP2=PI04 
END  IF 
TP3=SIN(TP2) 

TP4=C0S(TP2) 

TP1=TP3/LVP 

TP2=TP4/LUP 

TP22  »  TP1*TP1*TP2*TP2-2.0*CUVP*TP1*TP2 

TP1=TP3/LUP 

TP2«T?4/LVP 

TP23  =  TP1*TP1^TP2*TP2+2.0*CUVP*TP1*TP2 
iF(TP22.GT.TP23)THEN 

IF(TP23 . LT . 1 .0E-04*TP22)TP23=1 .0E-04*TP22 
LPP=1,0/SQRT(TP22) 
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LqP»1.0/SQRT(TP23) 

ELSE 

IF(TP22.LT. i .0E"04*TP23)TP22=1 .0E-04*TP23 
LPP=1.0/SQRT(TP23) 

LqP«1.0/SqRT(TP22) 

TP1=TP3 

TP3=TP4 

TP4=-TP1 

ENDIF 

DO  3516  >1.3 

3616  LPPV(J)=U(J)*TP4+V(J)*TP3 
C 

C  CALCULATE  CPPP 
C 

TP1=TP20/LVP 

TP2=TP21/LUP 

TP3» 1 . 0/SQRT ( TPl *TP 1 +TP2*TP2-2 . 0*CUVP*TP 1 ♦TP2) 

DO  3700  I=1.NLYRS 
TP4=ABS(LP*L0S(I) -TP3* (ZT-LOS(I) ) ) 

3700  ARG(I)»TP4*SqRT(TP21*TP21*UT(I)+TP20*TP20*VT(I) 

1  -TP20*TP21*UVT(I))/ZT 
TP1=CNNXMI(NLYRS.0.1) 

CPPP=- 1 . 0+CLIM/ (2 . 0*TP 1 *TP 1) 

IF (CPPP . LT . -0 . 9999) CPPP=-0 . 9999 
C 

C  CALCULATE  PERPENDICULAR  TO 
C 

T0PRP=1.0E-06 
DO  2000  >1,NLYRS 

2000  ARG(I)=SqRT(UV(I)*UV(I)*UT(I)*VV(I)*VV(I)*VT(I)-UV(I)* 
1  VV(I)*UVT(I)) 

TOPRP=CNN)lMI(NLYRS  .TOPRP) 

C 

C  CALCULATE  TWO  WAY  PATH  CORRELATION  COEFFICIENT 
C 

TWPCC=0.0 
DO  2002  I=1.NLYRS 
DO  2001  J=1.3 

2001  X(J)«VST(J,I)-LOS(I)*VRE(J)/ZT 
TP1»X(1)*U(1)*X(2)*U(2)*X(3)*U(3) 
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TP2«X(1)*V(1)+X(2)*V(2)+X(3)*V(3) 

TP3«6.67E-06*LOS(I)*SQRT(TP1*TP1*UT(I)+TP2*TP2* 

1  VT(I)-TP1*TP2*UVT(I)) 

2002  TWPCC»TWPCC+DP2(I)*CNNXM(TP3.I) 

TWPCC-EXPC-TWPCC) 

C 

C  INSURE  THAT  TWPCC  IS  LEGAL(>1 .OE-16) 

C 

IF(TWPCC .LT. 1 .0S-16)TWPCC«1 .OE-16 

CALCUUTE  OPT  AND  CQT 

DO  7100  I»1.NLYRS 
TP1=UV(I)*T0PRP-L0S(I)*LU/ZT 
TP2«W(I)*T0PRP 
7100  ARG(I)=SQRT(TP1*TP1*UT(I)+TP2*TP2*VT(I)-TP1*TP2* 
1  UVT(I)) 

TP7»CNNXMI(NLYRS.0.1) 

TP7-1 .0-CLli4/(2.0*TP7*TP7) 

DO  7200  I»1.NLYRS 
TP1*UV(I)*T0PRP 
TP2=VV(I)*TOPRP-LOS(I) *LV/ZT 
7200  ARG(I)*SQRT(TP1*TP1*UT(I)+TP2*TP2*VT(I)-TP1*TP2* 
1  UVT(I)) 

TP8«CNNXMI(NLYRS.0.1) 

TP0-1 .O-CLIM/ <2 . 0*TP8*TP8) 

CPT»LP* (TP7*TP21/LU+TP8*TP20/LV) 
CQT»Lq*(TP8*TP21/LV-TP7*TP30/LU) 

C 

C  INSURE  THAT  OPT  AND  CQT  ARE  L£GAL(<0.9998) 

C 

X2»CPT*CPT^CQT*CqT 
IF (X2.GE. 0.9998) THEN 
X2»SqRT(0 .9998/X2) 

CPT=CPT*X2 
cqT=cqT*X2 
EHD  IF 

c 

C  CALCUUTE  PARALLEL  TO  AHD  FINAL  TO 
C 
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TOPAR-ABS( (VTR(l) ♦LOSV(l) ♦VTR(2) *LOSV(2)+VTR(3) *L0SV(3) ) 
1  /L0SVA-VZ)/LPP**2+ABS((VRE(1)*L0SV(1)+VRE(2)*L0SV(2) 

1  ♦VRE(3) ♦LOS V (3) ) /LOSVA-VZ) /LP^^2 
T0PAR-6.6E-6^F/(T0PAR+1 .OE-20) 

IF (TOPAR . LT . TOPRP) THEK 
TO*TOPAR 

ELSE 

TO«TOPRP 
END  IF 
C 

C  CALCUUTE  X2.  THE  MEAN  SQUARE  LOG  AMPLITUDE  FLUCTUATION 
C 

X2=0.0 

XN(l)=-4.6 

XN(2)»-2.4 

XN(3)=-1.6 

XN(4)— 0.8 

XN(6)=0.0 

XN(6)*0.8 

XN(7)»1.6 

XN(8)»2.4 

DO  2160  I^l.NLYRS 

TP18«N(I) 

TP3«LO(I)^LO(I) 

TP1»TP3/LY2(I) 

TP2«(1.0-TP1)^^2 

TP1-(1.0*TP1)**2 

TP3-SqRT (TPl^O . 6*TP2) / (8 . 38E-06*F*TP3) 

TP2«TP2/TP1 

TP22»AL0G( ( 1 . 0^0 . 6*TP2) ♦♦2/(1 . 2^4 . 8^TP2) ) 

TP23*(LF(I)/L0(I))^^2 

TP24»NP(I)-TP18 

TP26*AL0G(L0S(I)^(ZT-L0S(I))^TP3/ZT) 

TP11»1.1-TP26 

TP10*0.6*TP22-TP26 

TP13»(0.69*TP22)/(TP10-TP11) 

L=16 

TP7— ALGG(TP23) 

XN(9)=TP7-2.4 

XN(10)»TP7-1.6 
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XH(ll)=TP7-0.8 

XN(12)*TP7 

XN(13)»TP7+0.8 

XN(14)»TP7+1.6 

XN(15)«TP7+2.4 

TP6»'2 . 0*ALOG (LI (I) /LO (I) ) 

DO  2153  J^.L 

2163  IF(XN(J).GT.TP6)GO  TO  2162 
J=L+1 
2162  L“J 

XH(L)»TP6 

TP26»2.0*TP26 

TP6=3.0 

TP8=AMIN1 (TPIO .TP7) -4 . 6 
DO  2168  >1.L 

2168  IF(XN(J).GT.TP8)G0  TO  2157 
2167  TP7»TP10 

TP14=EXP(TP8) 

TP3»TP5*TP8+TP26-TP18*AL0G(1 .0+TP14)-TP24*ALOG(1 .0+TP23 
1  *TP14) 

TP14=EXP(TP3) 

IFLG=-1 
TP17«0.0 
2160  TP1»XN(J) 

IF(TP1.LE.TP8)G0  TO  2161 
IF(TP1.LT.TP7)G0  TO  2162 
TP1»TP7 
J»J“1 

IF(IFLG)2171. 2172, 2172 

2171  TP26-TP22-TP13*TP10 
TP6*TP13^1.0 
IFLG*l 

TP7«TP11 
GO  TO  2162 

2172  TP26— 0.69 
TP6»l.O 
TP7*1.0E^10 

2162  TP16»EXP(TP1) 

TP16=TP26^TP6*TP1-TP18*AL0G(1 .O^TP16)-TP24*ALOG(1 .0^TP23 
1  •TP16) 
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TP12-EXP(TP16) 

TP16=TP16-TP3 

IF  ( ABS  (TP1»5) .  GT .  0 . 0001 )  TP14»  (TP12-TP14)  /TP16 
TP17»TP17*  (TP1-TP8)  ♦TPW 
TP14=»TP12 
TP8»TP1 
TP3*TP16 
2161  J-J+l 

IF(J.LB.L)G0  TO  2160 

2160  X2*X2+q2(I)*DP2(I)*(TP18-1.0)*TP17 
CALCOUTE  FO 

TP1=(1.0/LU**4+1.0/LV**4+2.0*(CUV/(LU*LV))**2 
1  )**(-0.2B) 

TP2*  < 1 . 0/LUP*  *4+ 1 . 0/LVP*  »4+2 . 0* (CUVP/ (LUP*LVP) ) *  *2 
1  )**(~0.25) 

TP7»TP1*TP2 
DO  2100  I«1,NLYRS 

TP1»SQRT(TP7*L0S(I)*(2T-L0S(I))*SQRT(UT(I)**2+VT(I)**2 
1  ♦0.6»UVT(I)**2))/ZT 

Q2(I)=DP2(I)*CNNX«(TP1.I)/(TP1*TP1*DL0S(I)) 

2100  ARG(I)=Q2(I)*UT(I) 

F0»F0IHTR(HLYRS . ARG . IDS .DIOS) 

DO  2120  I'-l.NLYRS 
2120  ARG(I)»Q2(X)iVT(I) 

FO*FO+FOIf}TR(NlYRS,ARG.LOS  ,DL0S) 

DO  2140  I==l.HLYRS 
2140  ARQ(I)=0.6*Q2<I)»0VT(I) 

F0-F0*2 . 0*F0I 8TR ( NLYRS . ARG . LOS . DIOS) 

IFCPa.GT.O.OTBEN 

FO**!  .48E-05*F*F/SqRT(F0) 

ELSE 

F0=10.0*F 
EHD  IF 

F0«F0/SQRT(l,O-Ci*Ci) 

C 

C  OUTPUT  LP,  LQ..LPP.  AKD  LQP  IH  METERS/SEC 
C 


LQ=TMPK*LQ 
LPP»TMPK*LPP 
LQP=TMPK*LqP 
IF (MODE. EQ.O) RETURN 
C 

C  CALCULATE  EFFECTIVE  SCALES  AND  INDICES 
C  USE  DOPPLER  WEIGHTING,  TP12=2.0 
C 

DP«TEC/(118.0*F) 

DL«DP/(6.3*F) 

DO  6741  I“1.4 
DLM(I)=0.0 
6741  DPM(I)=0.0 
TP12=2.0 
TP6=0.0 
LOEF=0.0 
LFEF=0.0 
LIEF=0.0 
NEF»0.0 
NPEF=0.0 
TP6*0.0 

DO  26CO  I=1,KLYRS 

TP4=SQRT(UV(I)*UV(I)+VV(I)*VV(I)) 

IF(TP4  .GT.  0.0) THEN 
R=»LF(I)/LO(I) 

RP»LI(I)/LF(I) 

QN=N(I) 

QNP=NP(I) 

TPS=SQRT(UV(I)*UV(I)*UT(I)+VV(I)*VV(I)*VT(I) 

1  -UV(I)*VV(I)*UVT(I)) 

TP20=DP2(I)*TP3**TP12*GWEIGT(TP12)/GWEIGT(0.0) 

TP6=TP6+TP20 

TP6=TP6+TP20*LOS(I) * (ZT-LOS(I) ) * (TP3/TP4) **2 
ARG(I)«TP3*L0S(I)/(TP4*ZT) 

LOEF=LOEF+TP20/TP3 

LFEF=LFEF+TP20*R/TP3 

LIEF=LIEF+TP20*R*RP/TP3 

NEf'*NEF+TP20*QN 

NPEF=NPEF+TP20*qNP 

ELSE 


ARG(I)  “  0.0 
END  IF 

2500  CONTINUE 

IF(TP5.LE.0.0)THEN 
KR=>>.0 
FR=0.0 
VEF=0.0 
RETURN 
END  IF 
NEF-NEF/TP6 
NPEF»NPEF/TP6 
TP6-TP6/TP6 

VEF»CNNXMI (NLYRS . XLEST) /TOPRP 

TP7=VEF/TP6 

L0EF=L0EF*TP7 

LFEF»LFEF*TP7 

LIEF«LIEF*TP7 

LCEF-AMAXKLFEF.LOEF) 

L1EF=AMIN1(LIEF.L0EF) 

R=LFEF/LOEF 

RP-LIEF/LFEF 

QN»NEF 

QNP»NPEF 

P2-TP6/  ( (VEF/LOEF)  ♦*TPl2*GlllEIGT(TP12)  /GWEIGT(0 . 0)  ) 

IF (MODE. EQ.l) RETURN 
C 

C  CALCULATE  RAYLEIGH  WAVENUMBER  AND  FREQUENCY 
C 

TP21«1 . 0* ( (0 . Oe^NPEF-O . 64) *NPEF40 . 38) *ALOG(LFEF/LOEF) 
TP21»1 .0/(1 .0-(HPEF-NEF)'*(TP21*LFEF/LOEF)**(NEF^NEF-2.0) 
1  /(NPEF-1.0)) 

TP20»TP2l*P2*(NEF-1.0) 

C 

C  APPROXIMATE  ESTIMATES  FOR  TP25  AND  VEF  WHEN  NOT  AVAILABLE 
C  DIRECTLY  FROM  THE  EFFECTIVE  SCALE.  INDEX.  AND  VELOCITY 
C  CALCULATION  ARE 

C  TP26=1 .46E^04*ZT*LP/(LOEF**2*LPP*F) 

C  VEF“LP/TO 

C  THE  STANDARD  VALUE  OF  THE  RAYLEIGH  CRITERIA  IS 
C  RC-0.1 
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c 

c 

c 

c 

c 


c 


c 

c 

c 

c 

c 

c 

c 


THE  PRASE  VARIANCE  DEFAULT  IS 
PND»0.026 

XN(I)  FOR  I  FROM  1  TO  8  IS  SET  ABOVE  AND  DOES  NOT  NEED  TO 
BE  RESET  HERE 

TP25*1 .46E+04*TP6/ (F*ZT) 

TP10=-0.49-ALOG(TP25) 

XN(l)*-4.6 

XN ( 1 ) =AMIN1 (TPIO . XN ( 1 ) ) -2 . 3 
TP17«TP20*TP26**2*EXP(3.0*XN(l))/3.0 
IF(TP17  GT.ROTHEN 

KR= (3 . 0*RC/ (TP20*TP25*  *2) ) ♦♦ (1 . 0/3 . 0) 

ELSE 

TP26=AL0G(TP26) 

TP18=NEF 

TP23=(LFEF/L0EF)**2 

TP24=NPEF-NEF 

XN(2)»-2.4 

XN(3)»-1.6 

XN(4)=-0.8 

XN(6)=0.0 

XH(6)=0.8 

XN(7)*1.6 

XN(8)=2.4 

L=1B 

TP7=-AL0G(TP23) 

XN(9)«TP7-2.4 

XN(10)=TP7-1.6 

XN(U)«TP7-0.8 

XH(12)=TP7 

XN(13)=TP7^0.8 

XN(14)=TP7*1.6 

XN(16)=TP7*2.4 

TP6*“2.O*AL0G(LIEF/L0EF) 


DO  41S3  >6.L 

4183  IF(XN(J).GT.TP6)G0  TO  4162 

J®L^1 


4162  L=J 


XH(L)=TP6 

TP11=1.1-TP26 
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TP13=-0.291/(TP10-TP11) 

TP25=2.0*TP26 

TP6=3.0 

TP8=XN(1) 

TP7=TP10 

TP14=EXP(TP8) 

TP3=TP5*TP8+TP26-TP18*AL0G(1 .0+TPU) -TP24*AL0G(1 .0+TP23 
1  ♦TP14) 

TP14=EXP(TP3) 

IFLG=-1 

J=2 

4160  TP1=XN(J) 

IF(TP1.LE.TP8)G0  TO  4161 
IF(TP1.LT.TP7)G0  TO  4162 
TP1=TP7 

J»J-1 

IF(IFLG)4176. 4172, 4172 
4176  TP26=-0,9808-TP13*TP10 

TP6=TP13*1.0 
IFLG=1 
TP7=TP11 
GO  TO  4162 
4172  TP26=-0.69 

TP6=1.0 
TP7=1.0E*10 

4162  TP16=EXP(TP1) 

TP16=TP26+TP6*TP1-TP18*AL0G(1 .O^TPie) -TP24*AL0G( 1 .0*TP23 
1  *TP16) 

TP16=TP16-TP3 

TP12*EXP(TP15) 

IF(ABS(TP16).GT.0.0001)TP14=(TP12-TP14)/TP16 
TP19=TP20*TP14* (TP1-TP8) 

IF((TP17^TP19).GE.RC)G0  TO  4173 

TP17=TP17^TP19 

TP14®TP12 

TP8=TP1 

TP3=TP16 

4161  J=J^1 
IF(J.LE.L)GO  TO  4160 
KR=1.0E*30 
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GO  TO  4171 

4173  TP16-TP16/(TP1-TP8) 

KR- (RC-TP 17) / (TP20*EXP (TP3) ) 
IF(ABS(TP16*KR) .GT. 0.0001) THEN 

KR»EXP (TP8) * ( 1 . 0+TP16*KR) ** ( 1 . 0/TP16) 

ELSE 

KR-EXP(TP8+KR) 

END  IF 
END  IF 


4171  KR-SqRT(KR)/LOEF 
C 

C  KR>RAYL£IGH  WAVE  NUMBER  FROM  FLUCTUATIONS 
C 

C  FIGURE  OUT  PHASE  NOISE  LIMITED  RAYLEIGH  VAVE  NUMBER 
C 

TP1-(LFEF/L0EF)**1.2 
IF(NEF.NE.NPEF)GO  TO  6105 
TP2»EXP(-6 . 4/ (NEF* (HEF+6 .4) ) ) 

GO  TO  6110 

6106  TP2»C(1.0+6.4/NPEF)/(1.0+6.4/NEF))**(1.0/(HPEF-HEF)) 
6110  TP1»TP1+(1.0-TP1)*TP2 
TP0«LFEF/SQRT(TP1) 

TP21»0 . 24* ( NEF- 1 . 0) *TP21 ♦ ( 1 , 0*6 . 4/NPEF) *P2*L0EF/TP1 • ♦ 
1  (NPEF-NEF) 

TP2*TP21 

TP3-1.O/L0EF 

TP4-1.0/TP9 

TP2-TP2/ (TP9«  *(2.0* (NPEF-NEF) ) •LOEF** (2 . 0*NEF- 1.0)) 

TP6»2.0-NPEF-NPEF 

TP6» (-PND*TP6/TP2) ♦* ( 1 .0/TP6) 

IF(TP6.GT.TP4)G0  TO  6227 
TP7— TP2*TP4*  *TP6/TP6 
TP2*TP2*TP9** (2.0* (NPEF-HEF) ) 

TP6«2,0-NEF-NEP 

TP6-(TP4**TP6* (TP7-PND) •TP5/TP2) •* ( 1 . 0/TP6) 
IF(TP6.GT.TP3)G0  TO  6227 
TP7»TP7*TP2* (TP4* ♦TP6-TP3* ♦TP6) /TP6 
TP2-TP2*L0EF» • (NEF*NEF- 1 .0) 
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TP6-TP3+ (TP7-PND) /TP2 
IF(TP6.LT.0.0)TP6-0.0 
6227  IF(TP6.LT.KR)KR*TP6 
FR»VEF*KR/6.3 
C 

C  KR-FINAL  RAYLEIGH  WAVE  NUMBER 
C  FR«RAYLEIGH  FREQUENCY 
C 
C 

C  DO  DYNAMICS  QUANTITIES,  THAT  IS.  CALCUUTE  THE  THREE  SIGMA 
C  VALUES  FOR  DOPPLER.  DOPPLER  RATE.  GROUP  DEUY.  ETC. 

C 

IF(KR.LE.O.O)GO  TO  6742 

TP26*HEF-0.5 

TP22«L0EF*L0EF 

TP24-NPEF-NEF 

TP23-LFEF*LFEF/TP1 

L»16 

TPS— ALOG(LOEF) 

XH(6)=TP3 

TP6*0 . 5*AL0G(TP1) -ALOG(LFEF) 

XN(12)-TP6 
DO  6770  J=l,3 
TP7-O.4*FL0AT(J) 

XH(J+S)»TP3*TP7 

XN(J^12)=TP6*TP7 

M»6-J 

XN(H)-TP3-TP7 
6770  XN(M^7)»TP6-TP7 
TP8*=-AL0G(LIEF) 

DO  6707  I-6.L 

6707  IF(XN(I).GT.TP8)G0  TO  6713 
I«L^1 
6713  L=I 

XN<L)-TP8 

TP3-AL0G(KR) 

XN (1 ) -AMINl (XN (2) . TPS) -4 .6 
DO  6260  M>1.4 
I-2*(M-1)*1 
TP7-XM(l) 
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TP6«EXP(2.0*TP7) 

TP6-TP7*FL0AT (I) -TP26*AL0G (1 .0+TP22*TP6) -YP24*AL0G( 

1  1.0+TP23*T?6) 

TP17»EXP(TP6) 

J=2 

6266  TPl-XH(J) 

IF(TP1.LE.TP7)G0  TO  6271 

IF(TP1.GE.TP3)TP1«TP3 

TP16«EXP(2.0*TP1) 

TPlE-TPl*FL0AT(I)-TP26*AL0G(l.O>TP22*tP16)-TP24*AL0G( 
1  1.0+TP23*TP16) 

TP16*EXP(TP16) 

TP12=TP16-TP6 

IF(ABS(TP12).GT.0.0001)TP17=(rP16-TP17)/TP12 
DPM<M)=DPM(M)-HP17*(TP1-TP7) 
IF(ABS((TPl-TP3)/TP3).LT.i.OE-3)GO  TO  6260 
TP17=TP16 
T?7»TP1 
TP6»TP16 
6271  JaJ*l 

IFCJ.LE.UGO  TO  6266 
6260  DPM(M)«TP21*DPM(M) 

DO  6306  M»1.4 
J»M-1 

6306  DPM(M)=VEF**J*SqRT(DPM(M))/2,i 
6742  DPM(1)»6.3*DPH(1) 

DLM(1)=DPM(1)/(6.3*F) 

DO  6330  J=2.4 
6330  DLM(J)»DPM(J)/F 

IF(M0DE.Eq.2)R£TUaH 

EHD 

FUHCTIOB  GWEIGTCM) 

C 

C  THIS  FDHCTIOH  IS  A  WEIGHTING  FUNCTION  FOR  CALCULATING 
C  EFFECTIVE  SCALES  AND  INDICES. 

C 

REAL  M 

COMMON  /CKDATA/qN.qNP.R.RP 
TP3*M^2.0-2.0*qN 
IF<ABS(TP3).LT.l.OE-4)GO  TO  10 
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TP1-(1.0-R**TP3)/TP3 
GO  TO  20 

10  TP1»-AL0G(R) 

20  TP2«2.0*qHP-M-2.0 

IF(ABS(TP2).LT.1.0E-4)GO  TO  30 
TP2» (1 . 0-RP**TP2) /TP2 
GO  TO  40 

30  TP2=-AL0G(RP) 

40  GWEIGT*! .0/ (M+1 .0) ♦(TP1+TP2) /R**TP3 

RETURN 
END 

FUNCTION  FOINTRCNLYRS.ARG.Z.DZ) 

C 

C  THIS  FUNCTION  INTEGRATES  THE  LOS  INTEGRALS  FOR  SOLVING  FOR  FO. 
C  THIS  IS  A  TRICKY  INTEGRAL.  DO  NOT  FOOL  WITH  THIS  ROUTINE. 

C 

DIMENSION  ARG(1).Z(1).DZ(1) 

F0INTR=0.0 

A8»0.0 

ZT»Z(NLYRS*1) 

DO  I  I=1.NLYRS 

A6»(Z(I)-0.6*DZ(I))/ZT 

A6-1.0-A6 

A7=DZ(I)/ZT 

FP»ARG(I) 

FOINTR-FOINTR^FP* ( ( (FP* ( ( ( A7/18 . 0*0 . 2*A6-A6/7 . 5) *A7*A6*A6 

1  /12 .O^O . 26*A5*A6-0 . 6*A5*A6) ♦A7^A6*A6*A6/3 .0-2 . 0*A6* AB^A6 

2  /3 .0)  ♦  A8/3 .0)  ♦A7^0 . 6*FP*  ( A6*A6)  ♦*2-A6* A8)  *k7*kB*kQ*m  *A7 

1  A8»A8*FP* ( (A7/3 .0^A6) * A7+A5*A5) •A7 

F0INTR*2.0*F0INTR*ZT**4 
RETURN 
END 

FUNCTION  CKNXMKNLYRS.TPX) 

THIS  IS  A  CRITICAL  RUCTION.  IT  ITERATES  OVER  A  LOS 
C  INTEGRATION  TO  FIND  DECORRELATION  DISTANCES.  TIMES.  BTC. 

C  DO  NOT.  I  REPEAT,  DO  NOT  HAKE  CRIT  GREATER  THAN  O.OOl. 

C  THIS  WOULD  COMPROMISE  THE  ACCURACY  SPECIFICATION  ON  THE 
C  ITERATION.  SMALLER  VALUES  DO  HOT  BUY  MUCH  EITHER. 

C 
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PARAMETER  (CRIT»0 . OOi . EM1=0 . 36787944 . EM2«0 . 63212056) 
PARAMETER  (NIP»il) 

DIMEMSION  AR6(NIP).DP2(NIP) 

COMMON  /CHNC0M/CLIM.ARG.DP2 

TP3«0.0 

TP4-0.0 

TPl-TPX 

DO  1  I-l.NLYRS 
IF(ARG(I).GT.O.O)THEH 
TP4-TP4+DP2(I) 

TP6«TP1*ARG(I) 

TP3»TP3+DP2(I) ♦CHNXM(TP6 . I) 

END  IF 

1  CONTINUE 
IF(TP4.GT.1.0E-4)GO  TO  10 
IF(TP4.LE.O.O)GO  TO  50 
CLIM>EM2*TP4 

GO  TO  20 

10  CLIM»'AL0G(EM1^EM2*EXP(-TP4)) 

20  CNNXMI*SQRT(TP3/CLIM)/TP1 

IF(CNNXMI.LT.1.0E-08)GO  TO  60 
CNNXMI=1 .0/CMNXMI 
4  TP4=0.0 

DO  2  I^l.HLYRS 
TP6-CNNXMI*ARG(I) 

2  TP4»TP4*DP2(I)*CNNXM(TP5.I) 

TP4-TP4/CLIM 

IF(ABS(TP4-1 .0) .LT.CRIT) RETURN 
TP6- ALOG (CNNXMI /TPl ) / ALOG ( TP4/TP3) 

IF(TP6  .GT.  1.6)TP6-1.6 

TPKNMXMI 

TP3-TP4 

CNNXMI*CNSXMI/TP4**TP6 
GO  TO  4 

50  CNNXHI-l.OE^a 

RETURN 
END 

FUNCTION  CNNXM(X.IFLG) 

FUNCTION  CNNXM.  VERSION  3.0.  4  APR  89 
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c 

C  THIS  IS  THE  PHASE  STRUCTURE  FUNCTION  FOR  THE  VITTWER-KILB  POWER 
C  SPECTRUM.  WHEN  USED  TO  CALCUUTE  THE  SIGNAL  DECORRELATION 
C  DISTANCE.  THE  MAXIMUM  ERROR  IN  THE  DISTANCE  IS  FIVE  PERCENT. 

C  THE  MAXIMUM  ERROR  IN  THE  STRUCTURE  FUNCTION.  ITSELF.  IS  ABOUT  TEN 
C  PERCENT. 

C 

C  THE  INPUT  VARIABLES  ARE: 

C 

C  LABELED  COMMON  BLOCK  /CNDATA/ 

C  (THESE  VARIABLES  ARE  USED  ONLY  WHEN  IFLG  <  0) 

C  3«N-2  =  INTERMEDIATE  SCALE  SPECTRAL  INDEX 
C  (N.GE.1.6.AND.N.LE.2.0) 

C  2«NP-2  =  SMALL  SCALE  SPECTRAL  INDEX 
C  (NP.GE.2.0.AND.NP.LE.4.0) 

C  R  »  RATIO  OF  FREEZING  TO  OUTER  SCALE 
C  (R.LE.l.O.AND.R.GE.i.OE-07) 

C  RP  »  RATIO  OF  INNER  TO  OUTER  SCALE 
C  (RP.LE.R) 

C 

C  FORMAL  ARGUMENTS 

C  X  •  DISPUCEMENT  DIVIDED  BY  OUTER  SCALE 
C  lABS(IFLG)  -  INDEX  OF  LOS  INTEGRATION  POINT 
C  (lABS(IFLG).LE. NIP. AND. IFLG. NE.O) 

C  IFLG  <  0.  CALCULATE  AND  STORE  INTERMEDIATE  RESULTS.  THIS 
C  IS  USED  WHENEVER  N.  NP.  R.  OR  RP  CHANGES  FOR 

C  SOME  VALUE  OF  lABS(IFLG). 

C  IFLG  >0.  USE  STORED  INTERMEDIATE  DATA.  ASSUMES  THAT 
C  N.  NP.  R.  AND  RP  HAVE  HOT  CHANGED  FOR  CURRENT 

C  VALUE  OF  lABSClFLG}  SINCE  LAST  CALL. 

C 

SAVE 

REAL  N.HP.L2 
C 

C  NIP-NUMBER  OF  LOS  INTEGRATION  POINTS 
C 

PARAMETER  (NIP-11  .E2—7. i609297E-01) 

PARAMETER  (D3*6 . 1788469E-01 . C3=9 . 5310179E-02) 

PARAMETER  (D4*l .3862943.C1—6.9314718E-01) 

PARAMETER  (E7-2.8490703E-01 ,C4-6.66S1693E-01) 
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DIMENSION  XPSD(12*NIP) .FP(NIP) .E4(NIP) ,JP(NIP) 
DIMENSION  E6(NIP) ,L2(NIP) .IA7(NIP) .JM(NIP) 

DIMENSION  XPSDP(12*NIP) .XPSDM(12*NIP) .XH(12*NIP) 

COMMON  /CNDATA/N.NP.R.RP 

I-IABS(IFLG) 

JSTART«(I-1)*12+1 

IF(IFLG.LT.O)THEN 

C 

C  DO  kZ^dO**2 
C 

A3=R**1.2 

IF(ABS(N-NP)  .LT.O.OODTHEN 
A6=EXP(-6.4/(N*(N+6.4))) 

ELSE 

A6=((1.0*6.4/NP)/(1.0+6.4/N))**(1.0/(NP-N)) 
END  IF 

A3=A3*(1.0-A3)*A6 

C 

C  INITIAL  NORMALIZATION  CONSTANT 
C 

A6=l . 00+ ( (0 . 06*NP-0 . 64 ) *NP^O . 38) ♦ALOG(R) 
A6*(A6*R)**(2.0*N-2.0) 

A6=l . 0/ ( 1 . 0- (NP-N) *k6/ (NP- I .0) ) 
FP(I)=0.24*A6*(N-1.0)*(1.0+6.4/NP)/A3**(NP-N) 

C 

C  CONSTANTS 
C 

L2(I)*R*R/A3 

C  XO=ALOG(1.0/SQRT(L2)) 

XO*-O.6*AL0G(L2(I)) 

IF(X0.LT.2.05)THEN 

A3®0.95i220*X0-l.0 

XN(JSTART*l)*A3-2.3 

XH(JSTART*2)®A3-1.6 

XN<JSTART^3)=»A3-1.0 

XN(JSTART*4)=»A3-0.6 

XN(JSTART^6)=A3 

XH(JSTART*6)=A3^0.6 

XH(JSTART^7)=A3*l.O 

XH(JSTART*8)»A3*1.56 
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JEND=JSTART+9 

XN(JEND)»A3+2.46 

ELSE 

XN(JSTART+1)=-1.36 
XN(JSTART+2)=-0.65 
XN(JSTART+3)=-0.05 
XN(JSTART+4)=0.46 
J=JSTART+6 
IF(X0.GT.2.86)THEN 
XN(J)-1.36 
XN(J+1)-X0-1. 49999 
XN(J+2)»X0-0.60 
JEND=J+6 

ELSE 

XN(J'H)=X0-0.6O 

XN ( J) =0 . 6* (XN ( J*l) ♦XHC J-1) ) 

JEND=J+4 
END  IF 

XN(JEND-2)=X0-0.10 
XH(JEND-l)=X0+0.46 
XN(JEND)=X0^1.3S 
END  IF 
A4=-AL0G(RP) 

807  IF(XN(JEN0).LT.A4)G0  TO  812 

JEND=JEND-1 
GO  TO  807 

812  JEND^JEND^l 

XN(JEND)=A4 
IA7(I)=JEND 
XPSDM(JEND)*0.0 
JMCD^^JEND 
E4(I)=N-0.6 
E5(I)^NP-N 
XPSDP(JSTART)“0.0 
JP(I)=JSTART 
XPSD(JSTART)®0.0 
JS=JSTART*1 
DO  900  J=JS.JEND 

A3=EXP(2.0*XH(J)) 

900  XPSD(J)=E4(I)*ALOG(l.OU3)^E5(I)*ALOG(i.O* 
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1 


L2(I)*A3)-XH(J) 


C 

C  2*SIH(K*X/2)**2  -  C1*(K*X)**2  ,  K<»DO 

C  «  C2*K**E2  .  D0<K<-D1 

C  »  C3  ,  DKK 

C 

C  Cl  -  COEFFICIENT  FOR  SMALL  (K*X)**2  «  0.5 
C  D3  «  PEAK  VALUE  OF  APPROXIMATE  FUNCTION  -  1.856 
C  DO  «  K  AT  APPROXIMATE  FUNCTION  PEAK  -  SqRT(03/(Cl«X**2)) 

C  D1  «  K  AT  END  OF  OVERSHOOT  «  D4/X.  D4  -  4.0 
C  C2  -  COEFFICIENT  FOR  INTERMEDIATE  K  -  D3/D0««E2 
C  C3  -  APPROXIMATE  FUNCTION  AT  LARGE  K  «  1.10 
C  E2  =  EXPONENT  FOR  INTERMEDIATE  K  »  ALOG(C3/D3)/ALOG(D1/DO) 
C 

C  D3»AL0G( 1.856) 

C  C3«ALOG(1.10) 

C  04«AL0G(4.O) 

C  E2»(C3-D3)/(D1-D0) 

C  E2»(C3-D3)/(D4-0.5*(D3-C1)) 

C  C1»ALOG(0.6) 

C  C4»0.5*(D3-C1) 

C  E7»E2^1.0 

C 

END  IF 
CNNXM*0.0 
IF(X.EQ.O.O)RETURN 
XX«ALQG(ABS(X)) 

X2»X*X 

C  D0«O.5*(D3-XX-XX>Cl) 

C  D1-D4-XX 

C  E2-(C3-D3)/(D1-D0) 

C  E7-1.0^E2 

D0-C4-XX 

XN( JSTART) »AMIN 1 (DO . XK ( JSTARt^ 1 ) ) -6 . 6 
XPSO( JSTART) »-XN (JSTART) 

J*IA7(I) 

IF{XN(J).GT.DO)THEH 

D1«D4-XX 

C2»D3-E2*D0 

IF(XN(J).GT.D1)TH£N 
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JE»J-1 

JT-JE+JSTART 
DO  800  JJaJSTART.JE 
J=JT-JJ 

800  IF(XH(J).LT.D1)G0  TO  850 

860  J»J-H 

IFO.LT.  jM(I))THEH 
JE*JM(I)-1 
JT«JE+J 

A4*C3“XPSD(JE*1) 

F0*EXP(A4) 

DO  870  JJ=J.JE 
JS*JT-JJ 
FN=FO 
A6=A4 

A4=C3-XPSD(JS) 

F0=EXP(A4) 

A0»A6-A4 

IF(ABS(AO) .GT.0.0001)FN«(FH-F0)/A0 
870  XPSDM(JS)aXPSD»l(JS+l)+FH*(XN(JS^l)-XH(JS)) 

END  IF 

CNNXM»XPSDM(J) 

CT»EXP(D1*DI) 

CT«E4(I)*ALCG(l.O*CT)*£B(I)*ALOG(l.O+L2(I)*CT) 

A6=C3^01-CT 

F0»EXP(A6) 

A4*C3-XPSD(J) 

A1>£XP(A4) 

A0®A4*A6 

IF(ABS(AO) .01.0. 0001 ) A I ^ (Al-FO) /AO 
CH8XM«CifNXH^Al  *  CX8  (  J)  -01 ) 

ELSE 

Dl»X8<J) 

CT*EXP(DI^D1) 

CT=E4(I)*AL0G(l.O^T)*ES(I)*AL0G(l.O*L2(I)*CT) 

A6*C24E7*DI-CT 

F0*EIP(A6) 

END  IF 

876  J»J-i 
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IFaN(J).LT.DO)GO  TO  876 
FN=FO 

A4=G2+E2*XN(J)-XPSD(J) 

F0=EXP(A4) 

A0=A6-A4 

IF(ABS(AO).GT.0.0001)FN=(FH-F0)/A0 

CNNXM=CN1JXM+FN*(D1-XN(J)) 

D1=XN(J) 

A6*A4 
GO  TO  376 

876  CT=EXP(DO+DO) 

CT=E4 ( I ) *  ALQG ( 1 . 0+CT) +E6 ( I) *ALOG ( 1 . 0+L2 ( I ) *CT) 

A4=C2+E7*D0-CT 

A0=A6-A4 

IF (AES (AO) . GT . 0 . 0001) F0= (FO-EXP (A4) ) /AO 
CI<N]CM=CllNXM+FO  KDl-DO) 

A6=C1+3.0*D0-CT 

A4=C1+2.0*XN(J)-XPSD(J) 

A1=EXP(A6) 

A0=A6-A4 

IF ( ABS ( AO) . GT . 0 . 0001 ) A 1= ( Al-EXP ( A4) ) /AO 
CNNXM=CKNXM+X2*A1*(D0-XN(J)) 

EHD  IF 

IF(J.GT.JP(I))THEN 

JS=JP(I) 

A6=C1+2.0*XN(JS)-X?SD(JS) 

FN=EXP(A6) 

JS=JS+1 

DO  970  JJ=JS,J 
A4=A6 
FO=FN 

A6=C1+2.0*XN(JJ)-XPSD(JJ) 

FN=*EXP(A6) 

A0«A6-A4 

IF ( ABS ( AO) . GT . 0 . 0001 ) F0» (FN-FO) /AO 
970  XPSDP(JJ)= ;PSDP(JJ-1)+F0*(XN(JJ)-XN(JJ-1)) 

JP(I)=J 
END  IF 

CNNXM=CNNXM+X2*XPSDP(J) 

CNNXM=FP(I)*CNNXM 
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RETURN 

END 
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